


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1990 


Thermal-hydraulic transient analysis of a 
packed particle bed reactor fuel element 


Casey, William Emerson 


http://ndl.handle.net/10945/28606 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 


Downloaded from NPS Archive: Calhoun 


Calhoun is the Naval Postgraduate School's public access digital repository for 
| (8 D U DLEY research materials and institutional publications created by the NPS community. 
«ist sia Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NY KNOX appointed — and published -- scholarly author. 

ia) LIBRARY Dudley Knox Library / Naval Postgraduate School 

411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 



























































































































avy stee ig atatte siti We 

















os i: ot JB OA A dite Rh wed Boe kbs Mm 
ea | m4 Le s a & ¥ F & a wey ond + “yi Wy eoik ‘ By } i 4 ae Pretty) ae 34, be A a ats 
I oot eee mete ’ a a a eT fy DF a a 4 Pa Vite ide tAt7A24. a My WSoi fal 
4 wi pe , 1 F ree r] A ey -_? x! 1 ® ais tye ¢ « a te 4 asst ‘ a? iy ' qf i ai “ or Are pti Oat ts S-a5¢ ot) bok 2 
"a ‘ : ier Ss ete” , seal : ns RtAnys Nerhpins reht. Ba ares ate Rt = aaah 
en oe | an in, iat eee ee tne S to} on es 2a settee Basie bis Se th fer oe 
i Poe ee gato set's be - ed steal Rn Ks Mita + “ i, Asyyp ate | we aiceely. “a on Bea fates 
were a Sarre Rac gt a eee Ra la ; 3 te dae 4° eh Biaea uae 
- ; > Nerina: os : 5" ‘ ~ ae Sige i “ 4 vit : vr t ify * 4 o s% Px ord. Kid ne a ice 4 r ar Catacd sat Fri a ii ee 
EB C F LO i F o4 Ms 8 PA ae 8 2S tae a vay ey i Alase's “ack, 
ti Je 7) A s eae fags Meh, ae See fee E¥ ares ge to We ra iy Keil fiag iis s eee a aes 
Pa? a ; eS 4 a “ae "\ Is iv ans ef 4 : ‘ vy fs r rae ope. A . ie ua vhy i “og ae , He q es iavee Ane eer 
e “a ‘ 5 nH ts AU what via A ots 7, OI Fade G § ies LP tk. 4 
ae i ho @h8 4) \ RE | - wee HA ; o¢a 2 vA . a rok fs 14 “ 7 i] 4p 
' ete 3° ¥e te " t * % 1 Conep oa be a ee aes ‘ 
5 + bd s . ‘ Le és Ba r 
: 3 ' t Py 4 Phy 6 A Le * an Md iF" "e y 
is 0 ue +4 ‘ 1. I pier say 8 R ; yt ‘} y ahe ye a be 
i - = ' | » jay ef 3 vs te z 
, be ib ike are té ' ae e Yi gun: ‘ hd ye Ry u,. 
i ' ar 4 i ty ead) oe ee ae a ve Ne ye fat Loree, F 
' ' ' : ' a ts 4 hae ; 5 ‘ , rire A aye or iP is 
. o ie y . eee, he ' ’.,! A “s ats a anaee fi Sedu a ay Me vat y fe Ay ished pone a Te : ; fe Deets st 
ay ‘ " e "6 ha rik > j oF . e146 = cara: # % x . le 3" 4 | *2 Ni. hs  Faagied: ‘s ” fatitaga areata hss : ET pint ps 
ae U : i oft tae de to & 4% be 3 Ftd c yy . , ‘ng gt rT Fiy Lae Aan eh aM MEN Rh wth 
: ar te * ' Veg gen oth 8 Hh : $ “? Bt rer hse oak + Sat4s 
' ; ee ke ee Me wy se PSE alge teeg at af eity Veneta aes pee 
. . a a oe a ry ‘ a =! hy bee va bee ‘eens’ § chute) 3 ae wi ey ays ed's «fe cy 
Cie ' ao =, Fer Jay 9 is = etad are i avid all ls my -4 wre ‘ i val 4a Se, aK rhea ath ot Bite ety 4 ethene ott Sea at aoe ON et Peg 
; ; A pa ™% wee aire Wie LY 5 i ; aoe ’ Ve 8,8 pth era la ’ ppl Ey . ui rt ta, yy = 
ee ga ur 7 ual , a " 4° ei, woe : ae . 4 j Nt ae mn ae hh ‘ a af. “ pane ‘ Wi ee vat rn A Mi oe ee eo a ae PaxY ee i 
y oe ry us ' rate cera at Ve ‘ . t x oe : {> ; O v, ' ty tate . bg a¢ ni? “9° ay é t: ast AO ip in td “ee bey “aby sae) % (area ee a ey cl f 
5 Me th GED ‘ ‘ oe Bory ore a : f i an . Petes 4 4. at vy agi ue { ae owe re rae 9 iii aha a nath des ai tig hye ea s pay ze 
rear at ee ‘ ee 6 st agat oO seeders Pr eh ee A # A a ( AN a sai? H head. 
ae U ‘¢ ‘ ‘ 7% 5 ‘ee Ak o s} We iyi" $ as “4 ct tee ist od Sh, ww? f, Wy ‘a! > daeetiis 
ve ay ; ate oe . . ; Cys ae eT AT i Lary be tr een zeae pss) 1% ti sp ‘Ag pak z. ts Gh Aina bine iB te art a a Ot FD ores e bes oe y a pie, Pipe 
cy, > i eee ‘ Sa La ve ee CAMO ¢ 2A tee «flaatle es ide is td EM Teg wr) ‘Tie oat sae eee ee 
ee ps pea Bie ante es wie th tegdawe + ONIN aay! Sass ie eit arta seeps ote eye aa net is bor Sa Set ss eve poe sa: 
miased peg 1  enann m4 ba Tt ete + ote . ’ 4 $ rf Sy ys s iA Macau ‘i 0 i 7 sae ety 
eae is : abet ald ate Sy we her nee i" fagoin 4:47 2a hore Ase : 
: : i : : 0 68 TWh go LA . ie aes Ca ¢ 4 : ae 1, ne ‘wt tie * mE “a : ru ag ras sh he & Fs qi ss id ad) ert Ae: i haha f daar a ee ae i cot eh at . srs 
' 2) ue ene renee let 2 ae’ tat yess ahivitcsd) <0, % wae Ren Ae. ; oh oR Heed ~<A Me Wy ay a “A eas Sy feisava Wreck are heagetacu se 
aie a pees hee Oy ee Re ’ eel ae i mee ovata whi, 4, tal 4. 3 Dy jute wg a aye vhs Tali iF 4 bs, asd aking fe hats : vied es tee “aig ips iege sips. ay alas vont sy. Ls BARS PEE RIN gs oA fd 
mene at aL Mee ert oe ei tad Grint ed, aceat te is Heit yh wees bce hate itty Ae af z 1m 586. 5 yeaa iP Tae ra hire ees 
bee 8 “4 ' 4 Ai ye oo A PEPE Usié eek wre rs ten thane seis Hain be ai Hi saci aOT & 4 i54 ee i ? 3 é BA caie: 
v4 Of ve 14 aT ‘ e STC) F + {,, Arua : ite 4] “43 <3 a's 9! & taf ae 4 ‘ “ina fot Nisan hg re) te hl fF) sae RS, Aiped aang Pots ier tas eae 7) feta 
Sees Migirans Al ery sy cit tie a se gtpats “Eh ee, | iy “pect o2 hasaca $s yee Ay dassiatt An fede: So r ei 2 (lesa ae (4M 2 eit oan 
Dot: er Be OR tens he ost Gk g ae 7 “y's rea Cay ae wee aceite ASRN gh SUF ‘alate ar es right Seas Phy chads avin, e weak 
i Rei ae ‘ ie ef {Pau re ' bes ate oe ts ¥: fel sit 4 pa TH ih a oN Fake ” Boy “iis are bes eee ae eat ies - ZH 
F ' ; ner mes ae hah ‘She neta oS fh y 93 4H ay hada agis CLA s 
: Ab te ae! TY : Rae ¢ vie i 5 yA’ oe A uy’ aca hg ‘> oe Thy) sed we Ass +) be. “he att ate? CUES sth hen Ares bee yoy afc ng ts wh ny oF es cee Sealy 4k it soarialee 
1 . s H 9 A’ é ' £ 5 s, “ ote 1‘ 2 AoA ’: < aa 5 wiry, ¥ rae 
Bee ee van eee] 7. te : ; +: Eg Uatiane os ta a ae rite onl 5 A of Se ie al ag pie ry yet ys ny, oe ai woes haat igdine thgadueusee a in i aiein aha" ey ie sigpak ree i 
: ] ag be A ' LP whe 98 ° 2 te! ye eth ee ny wank £ mot i t": ‘s *}s a. re thal (aad, 4 ats’ ' 7 a) a eho * MR: fie Vth ¢ it 2 aapettt tetra Bote er a 
ics Lee uy ‘ e ae® orb apace @ ORE eae ew nett ee : ie 
a . Ys 2 oe ‘ Pa Hes oh tha?) ob aslen At thee {' int ax sie “ad J: a “ay wr retain 
i ny Te tart sa oy at ¥ ees asl be Tatas Ars eae Ser 
ceeguipis 8 ary rit tee pate rudy” nat ty CO $i. wh silly 4 veut 29s diate oe fergts anes 3M x: ok Maas cy 
ee eft PRD Day atk 1 rath We Bost f 
we t ‘ P qa vatalgs a By a ay fi a tee vena laid Ait Taig Wide oe Goh fee's, 4%. 
: 5h a ‘ ; i'l ; f aay Ye : dinty: cote af ' De Policies ite. “a i Chl ass Mie bee at 
Pf] chee Tor Cet ae ce 2 * iy 4 aa fate ay i STL ot tae Mh y baci fade. Par ty Ee ear SORTER tie . raah th vy ahs ip snsa or rr meade eben “an aeieg A Singh 4} 
Peau i eS er ae Le AL! ee ah aif toy Wig die sSyycatin paint aaa seit aS ba Chita ta wake, ig Fa aks ayes een athe wate 
‘ er] rs fuPefs aw 8 nati i" . ‘oy : Me Xe “ Mf j 
' Dette 9 Be here oat ay é Ys! 0h og A ut ag! ° Thaw. ashi tatty es ORS Rett bre Karaja no A a Hes siege Dyan), ih pats Set tae nt ae (ay. el Ne 
nee pues erent a ; say ie baad ede HS ey Ol anit % fot athe} ays ‘e ae ‘aie xs aay eae i 5g “bc 
: Pee ; yon fees Rt nits Serb ics at 1¢ 




















































































Sthrata; rss btwn tralh wpasie Pai cd 
i panes acalhy Tan ed one WEA oa GaP Cong GA cera paces er hgite 8 pedal : eee Bases 
Nee? 3% iQ ts ee r i “Uy"s ¥ ‘ae Le ati one te Oa Myiate cAih Be is Ee: yoda heal iN es wiry sara faye hel Deby eats ia ah EZ Tide 
: ee Abst 94, 8 tasauh wt adant Acaéh ia tea sea oe Tawa tg are cate ree cS; Ee RiSiete toad ands 
; ee obhetacw 54 ; 943 nae Piet | geAtd a Eby 2 Sa mali ork tires vias “tt ve ey aay Ha bere Orie arp Ht py hes 
| a idl SICAL Nae Ale, SP ub ates a ares bed 
: agse Siyierea ! ‘ Bat pact Pe ASI ' ° t ‘thd mS FACS AR ies 
aan eeeel dy ud auieein sal Ri oa 8 “h feted ata Ki ace t by , nies wis PU ARs vei At one 
' aa é H sel tare 2, *M At ge, 4 ; 4 5 ert) ' ! ¢ ; 
od ‘ ot ie Nee fia tage ‘ttele dear ie ee big! Wh mae (3) Peg wienan ne? ye he Salam bel hed Asean 
Cae Ss oy ar Oain rae inate Ce fu, gatathty pete fe: we mies vis eee ASME Id th Heelies ds Paes $ i / gg Bo lei mal vita gy 
‘ ‘ 2 . i ate rob O Tae i Rc y rH wy Pies Diy? Mees apmen 7a re] LA > a! td eur Bec Alag fore Seite ed : ey be 
ar UP de ai 53 ‘ oom bs “a as ie Leas Ne : : i ra ake :: often! we Te ate Avy: “ fe et ¥ HM tye wie oke ay Jas aeteis $6582 TESA fe eaitiie seu rice ‘ 
1 ae ee ea Md he ce x teen bin MY ee ; hy itty ee recas 384 fai ae we hig heats Ht Are 
: a PRL i ife's r fie fegr did ok’ hte ks 3 Me} 
' ’ . for wit 4 sa: ste ; Gs pt lc < 


‘ * ne 4! « 4 

reife ay pd a Pie “ nh Ae f 

a ate Bibeat RS en Sey, hh sees ie 
ath a ae Lae in ON Rake Sig, Sei NGA 

tas ab Ata tna Toe at 4G cae! % 

Mhdags heats, vA vase chl va \ebtebes pia dies tau de 
A th af 1 VNIE eet : ea ities HORE TT 

~ lof seit sen) {i hy were fh ie figO geht owes Heh age Fe 4 Oyen eG hy 5% Ubi fe 
| MA sgp ate fe rime k ts sy tes 45} Pe aet byl his Falnt gt és oe ategae ig At Mees jie A va 
“si ted oon hie Tatheuts Pre pra tS Hf ee , eS oe 

34 


7, 
Pt Cer ere Cota bate, gins Sq getiacats tiga ists ; ie 
: naa arta 
SS 6 ae ei 
Bi he nth 


rains 
oo 


iat Sr] ; 
, Maral, a Sa, st AA-32 itt 
a ate sia rng ie iN 
™ ‘ie og KE a Sis 7 eek af! Bia se 
VJ 
















it ‘s e oa 
























” 
ei 4 a fiat rts 4, Pad ite! aia 
: i » tet “fi yes 
ag? “ie os one ts tio 






= . ‘ efi leah 

AU PAC Pay re Ainge e aga Aik hate Cee 
' tub « Pers f iP ade pv 

Nie Pa DAL ate ¢e) 


385) 79/5 
‘ ae? 

























































































































ah m=, fh. » 4%e aid 
A 4 a’ afnne la i a (ait Stans oe WDA ty ey i eae bbe) a i 
1 tub HA dass dats ab ‘1 
iP nth thet eta ttie be RDU I ‘ ae ne a se 
a2 ett she A 1% Hiake ou rie Ava RENN ee eral pea Ye p us9 : £5, 
: Ac o 22 tintg ate we S ate oe ee oh pean cul = 
, rf t Ph¥ nm 7 % 2 eh tile 
aut rd te paras te styh wand eta Ae a + tek L cs ke 
; RAS Gay Rk eae) ign A 1584 ae 1 Pio a1 : 
pee, F850 Sin 5 aes ates mee tier ws AS tere 
re . ° af f SPs 2h Ei 5 th Mot ft Lice, | By 
cba) d 1 5 5 J f op Thy}! 
1 ay ott hk «an wa: ', ie ne he erate oat ROY ern re aime 43 
: ‘ an Cea fee ar ie oS beet? RN AR ee a (7% i ; “4 
: . aay Pree eae . * ’ 4 a°q0 » No rie Pr ie Efe MF saps laity has hk Pa ie 
' tae Fe s , ®t 8 ' es fn, { fe Fas a? Share fe Ont oh erie, tat 
so. : Las aeagne iat eget eee 
ee “ he pen eras ey ‘gale ake is aie Ssartniles 
ae eax ' oon ' sine ba vy 3 SiMe ‘ ee eas ee ae Meat sii fan cS tye Me; ies 
' pac oka ee ea - rng pee tft AS 4 Se eae seed 4 f 2 (3p he {4 
ab 14%; 
— 5 4} , pve ahi it nee 
; 2h a Phe, sts en tts ae 4 Aes 
e 443 f Bafay tite es aati tosrd Hr ies opp ; 
‘ Ha Af. 









43 z, 


aacids ares epintt 
Ln ; 


2 < 
bel BATTS ws 
ie Ey 23 jist es 
SNe jays: ast co bee +f Bho 
Srihari ets 7 eae mh ei i Per ct 
Baty | 






ayhd sisal: 
¢ sete seadagh 
nh Taratgta oY 

eee iy aide $47 5 













Nee ye At Sift at es ered a ae es 
Res ees ‘ iy anon ‘. Doe ate os - 














erry 


we (4ah | % tie Saray, 
ie i SH pe hee 


£444; 
Pe sat #.. ee? 
iJ 


<2 Sree 


wen 
re 
ite. 


ies oe eons 


AP 




















: cpghtiedl f ae “f it f oe is aie 






















mieite ie nae ‘i tit 
; e Ms Seca Re fie oe 
v ¢ hea are eet oo ft = 
= res Ane ‘ 4 pas fe Ag oth rik yang 
‘ Teles, stout at # a 
ei aieeteas MCE 2a ase 






































; vy ae an ae ae 
a 
boat 
1 S ae ae it ue ay 1s 
F H shy y tel Ses th aay bd 
, = Fe v3 Bits loan iS 158 ai é sp iS xs eos! AM IEy ss 
made porone Haute Hiatgl iat ee eek a Ao cox 
Gey) agiselem vy ale ts i H jait 0 tr) r af Nee asesd tne ies A$ ass 
neon Wy «4 x! Wea arr agents ny et atone soe tae i Sh ape crear oR Sra 
St ae i YON Sta oaaa tte alo ae ih se eats apie er; aa nex, ae eed 
th fe ae Re Sth UE at! lee aa 1S aS Rie a ont an aac a ate 
ee as , yeaa A «, RIGT ' fhe TE Fg ge pee poy 7 acu p s ’ o? pr 
ra ialaxit sanady, ot B Pyare He Tt (pitas Ye Neate WAN ap hat parrag en ban aedens oe tek 
P AS Gaey oe, Ret odes eRe als ieee os he ut SESE MAR AY 5 ESERE SEB 
os ‘hae paar si oped atte Feed as i ik ee “ eds Sh. tens done seke wap te i 
Ce POL GiNitre arte ne ae ete 
ast 175 nts e + 









ir f, 29 i593 4 

ERE ES | aptes * eat Re AG 
Bf ‘giy fie wy eae 

Re rte 


*) 
H "J . . Pika oy oe 
nae ase “Wyte Leash piled ifpala itt Nake fe 
ee Ai "yas sisiae is Sieg HS Ra. ta Mies ¥ 








Ms ar) ote ey we ‘4 ia te gts; ue 
an) fener Res vit et eee is atest APE: M4 BAL ve 









































































































pee 
a 2 ieee ee et 
c Sos st ite tan es iho TPE Seg oF eb yet 
% abiigiaceng “ at TE En Pat yeas He Ne : BE bt : 
Me ets Mo ee tag) 4 fit ye ty rity Hh sede H ae a ye y y 4 
: ty “¢ yee fi ee gated au ; f ai 
nt, ny rd ?" ae ate Use is ile a cas oiee ati ipa eee iets Ke pri ere 
tact digi i BOAR tid been G a ae ele ft Hs gtr rithan ey 
a ro Re ix! “it nae: i Arana pha Shs putea Ate sie pare gatae 
. Z ing 4 * air. a i+ yaoi | ‘ +, ry : ae ft} ms eT raid) Bent ¥ Fe. e . 
s ‘ sat Ant eat tauN ie id? 3 Te sy abe: bgt anya eae iti tee $Fie 
‘ ~) om, ates Sey fs 2 ERAS i eH a ent ret : bags fo 4qney aha iiss ig ee: oe a naa Wh Seer te ree a o2) Sop rate 
Thy one ¢¢ Ny ats, fe & aa © 0% “5 bd oct wise Fel a ts ea 0 MAG MG pe site ae seek :. UA ad ti ah Uitegoteg sys este. 
pee eine Se io oy ay reed rit th Chin's ing iA etsy aot iii et ; iit age une as hee Lasik om nee 
Sere pe? ae nat, and ectiza bf ie yas sed Rehist gi ; Bit : ait es a7 id ¥e atts sai 35 ae 
git ry la Mental be ‘ 4 hits, tyeem, ayedtyt pnt stent feats? ona ie my it es Peis shah’ fhe co geed pe 
ogi woe a eta an ee t f y ‘te, fAagitil, es Suey AUP AIRES oy Cmte at eae uategua's 8 wi epee 
Sy ANE haa i MORIN ER Me ae eet tm Gi A Re a al ee 
, a { . 2 Sete e eh ttt del eug Fy” 9 t ir en F ue 
eRe UME ay Ce sleca careataatiarigiltiole ay ee Anal Mane! 
= Oe eh t a! ? ai ute « ‘ ei ony Rae re tk ae ane Br G by J. sry a 
: east pie A oe Tied oleae alai gre ALPHA ICTS Mg CHAT UF wfset i 
to ay ; ih ao Be Aen . of Pas « OR gael ¥ Hie las ty ator Sithwal. 3 iG ena 


my Pr : “ 
ayia am dei Phi a oe as 1 Teed eeu WSO 


4 ae taste SEE @ . ny i det ten’ 7 


a Se Ae es 
Bak Bika Big: Senger Pe ete ane arenes aban 
te: oat eae wae Ne oa “is ane bd ace et Jt 

at perianal eal tats ain eran epee ey ( pe a ibaa brett 

"3 44 oy Ae 2) peat 


¥ 33 ate the f Ged eimeaes Prtetrtls ae Se heden 
Ae fe Ree a ras te 1s 4 fetes te laghista tre peat eat 
riety gt SAN RP ¥ ora LA (hese Ware On tae wee as . aroet 


* ad et 2 084 
Pet ab Santee 2 gaat ue 
otegs sgh moda! eve eae fay Hay 9818 a eS AAS, ey, eat age te 
wks Catt Pee PY asl Hades eT th iat ataha Sed gs ee 


: ernie has) 1 
Bho 6 vA puioe ya rat eboney 
one ay Jy  atiet. Todi hse Sener fiken Bat <3 Ree eas poate rubies Bs She tae tm 
; rm at Se Pasties Coaheal ts a: + 1ea He a hegiyar Mine We por! ms ab ee sta pee 
peels fase oat S ita ters ee 2 ae ghey a3 : + ¥ BAY ears ahs Neve. ty Le » Bhetrs ere ie Sioa ry; ts re 4 
+ A 71 pry Pa * Kat a * re fl . TERE 235 A 

‘ale (bast she “Aa wi ie syise 7) Pv moi gebiaies uewu gaye net : races : : - 

te 


: » ben 0 5, “ sat Sad 
C he ea Omer cae 














caseni tat driegs 


















teauteiy SOF Leh a ed 
oA aah apes ROE Wi rhe he ketnted nega 
** ahi mo oth Shiga a at yale ate > EG D\ Cia x Ot Pa : 
oe Me x arty as ah Set aft fyi Ae xf are ght APOE CE autos 2 ade 3 
oa 1 . we : be st. mh fl Tete fe nS. re x97 py Sto 4h Shs 
| ek ee ee! 
‘ ae a ee . ere: Be = os, Dit aes ieee ae eee ha ee Bit aetfl y PER phe ety! “im Che 
a 4 a . 4 H | arate ey 5 ay Wp: 9 Yo mig Mee? “é init PF es aed 
‘ ; ’ , or : oi *, Be ’ one Unsere stn ie oe u eet oy hers 
". oe d Us wee , Ess © SATE Fi os vane 1 tS AL PP 
oe bd £ . “4 
Rak i 


ais fae { BOB IRA! 13 
ane uf oi Wale, erates a 6 NPS SE B Se epee ts at, wt ae " afion: 
y : ats p a nly sue id Semis anny aie ae paisogt 3 ‘ to 
. > : , . : ‘hot ; ai : “t wry ei mae f, iets ie eat 4! 
é 48 "1 : aria | : r » F e 4 - > iy Vi thats, i 


one : 
a 















i GaSe An aL Ey NEA aera oa Pr tig 2 
ne ms hte a “PENAL Tite Alitere Zing: Jagive fe, digi SLE bith f tts fet A Le Neue sieite 
r) 4 é . a. = : : Ly S3z x 
> ’ ‘ . =! r : « Ph fap “¥° F r 7, h ie * 4A ". 
f E x : ' : ' , ¢ “a! 7) f of. ane re ; PUL ANTES aero diabyy Gy Ma: as See 4 ys : SSSeapa 
7 : A ia rt Meee Fae ; Ve ise omit valet Gectigt Al i th Pernee f 
roi SO pias toes ey meat sey : aa 40 eee 
* ' ® ' Sah) Sa eared At bh he 4c arnt 
are oe ie os ; 3 or Vols ey puch fs a es ey ach tar: errcid te ae 
; a. yond a bipry Beer SEN Ge ey ak MRR Coa oes 
eae 7 mn Peers i. E as a we Rati : eee gs Sah pie ae ees 
‘ e . Tk BRG uty TF Cat ews hia ate 
. are a ok: m4 He Ra cb iveth aero fal Ae thie oe Shar) 





sacar ee 
sae 
pn, = 

















4 e 
ed 


THERMAL-HYDRAULIC TRANSIENT ANALYSIS OF A 
PACKED PARTICLE BED REACTOR FUEL ELEMENT 


by . 
WILLIAM EMERSON ASEY 
B.S., Engineering Physics 
U.S. Naval Academy, Annapolis, MD 
(1976) 


SUBMITTED TO THE DEPARTMENT OF OCEAN ENGINEERING AND 
THE DEPARTMENT OF NUCLEAR ENGINEERING IN PARTIAL 
FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREES OF 


NAVAL ENGINEER 
and 
MASTER OF SCIENCE IN NUCLEAR ENGINEERING 


at the 
MASSACHUSETTS INSTITUTE OF TECHNOLOGY 
May, 1990 


© William Emerson Casey, 1990. All rights reserved 
The author hereby grants to M.I.T. And to the U.S. Government permission to 
reproduce and distribute copies of this thesis in whole or in part. 





THERMAL-HYDRAULIC TRANSIENT ANALYSIS OF A 
PACKED PARTICLE BED REACTOR FUEL ELEMENT 


by 
WILLIAM EMERSON CASEY 


Submitted to the Department of Ocean Engineering and the 
Department of Nuclear Engineering on May 11, 1990, in 
partial fulfillment of the requirements for the degrees of 

Naval Engineer and Master of Science in Nuclear Engineering. 


ABSTRACT 


A model which describes the thermal-hydraulic behavior of a packed particle bed 
reactor fuel element is developed and compared to a reference standard (Tuddenham, 
1989). The model represents a step toward a thermal-hydraulic module for a real-time, 
autonomous reactor power controller. 


The general configuration of the fuel element is similar in construction to a design 
studied by Brookhaven National Laboratory and Sandia National Laboratory. A bed of 
small (diameter = 500 um) fuel particles are packed between concentrically mounted 
retention cylinders referred to as frits. The element is cooled by parahydrogen which 
flows axially through the inlet and outlet plenums and radially inward through the fuel 
particle bed. 


The momentum integral approach used in the MINET code (Van Tuyle, et Al, 
1984) is applied to this model to balance the fundamental mass, energy and momentum 
conservation relationships. The element is divided into only three control volumes: the 
inlet plenum and cold frit define the first control volume, the fuel particle bed defines a 
second control volume, and the outlet plenum and hot frit define the third control volume. 
The solid phase of the particle bed is represented by a single node. 


This simple model was validated against the reference standard and compared 
favorably. As a demonstration of the model’s flexibility, a number of variations were 
analyzed. These included variations in fuel element geometry and the initial and final 
values of inlet temperature, inlet pressure, and outlet pressure. As a final demonstration, 
a cluster of nineteen, 1 meter long fuel elements, arranged to form a core, were analyzed 
for an up-power transient from 0 MWt to approximately 18 MWt. 


The simple model significantly decreases the time necessary to perform a single 
analysis. A transient of 10 s with a timestep of 10 ms, for example, takes approximately 
45 s of computation on a desktop computer equipped with an 80386 microprocessor. 
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CHAPTER 1 
INTRODUCTION 


1.1 STATEMENT OF PROBLEM 


A packed particle bed reactor (PBR) is a gas-cooled reactor similar in nature to a 
high temperature gas-cooled reactor (HTGR). Unlike the HTGR, however, the PBR 
packs small fuel particles between inner and outer retention elements, designated as frits. 
The PBR is appropriate for a number of gas-cooled reactor applications and, in particular, 
seems to be most appropriate for use in space because of its compactness, high outlet 
temperature, and wide range of delivered power. 

The PBR is proposed for use as a power source for a variety of systems in space. 
The requirements of these systems range from relatively low power levels, on the order 
of kilowatts, to multi-megawatt applications (A-1). These systems include nuclear pro- 
pulsion rockets as well as systems which may require large, rapid power transients. The 
proposed PBR’s will initially be unmanned and will therefore require an integrated 
autonomous control system. 

The Advanced Controls Group at MIT has worked with the Sandia National Labo- 
ratory to develop a control system for a PBR which is to operate in space. The system 
which is currently being considered uses a control algorithm based upon a real-time 
model which must accurately represent the reactor so that proper control signals are gen- 
erated. As a step toward the development of the real-tume model, a thermal-hydraulic 
model was formulated by Tuddenham (T-1) to calculate the behavior of a PBR fuel 
element during transient operations. This model was developed to provide detailed infor- 


mation of fuel element behavior and to serve as the standard for future models. The 





objective of this investigation is to continue the development of the model-based 
algorithm by evolving the thermal-hydraulic model of a PBR fuel element and begin cou- 
pling the element with other system components. 

The automatic reactor controller will rely on the thermal-hydraulic model for per- 
formance predictions based on temperature and pressure inputs from the system. This, in 
turn, will allow the controller to calculate projected fuel temperatures and allow operation 
without exceeding safety limits within the fuel. Further, the model provides information 
on the temperature and pressure of the hydrogen coolant which 1s used to obtain the 


amount of reactivity feedback from changes in coolant conditions during a transient. 


1.2 APPROACH 


The major emphasis of this investigation is the simplification of the thermal- 
hydraulic reference model provided by Tuddenham (T-1). Although his model produces 
a very detailed analysis, exceedingly long computational times are required to reach a 
solution. The simplified model allows for a method of analysis which is sufficient in 
detail to give an accurately calculation of reactivity feedback within the fuel and yet is 
sufficiently simple so as not to require excessive computation. 

Once the less complicated model is developed, it must be validated by comparing 
the model results with the reference standard. The reference standard provides calculated 
State variables for some steady state and transient conditions. 

In addition to building a simplified model of the fuel element behavior, preliminary 


steps are taken to couple the fuel element with other fuel elements within the reactor as 


—_ ss 





well as with other core and power plant components. Future developments can then pro- 
ceed with exact modeling of a complete system once specifics become available for each 


separate component. 


1.3 CONCERNS 


The major concerns involved with the simplification of the reference standard are 
those which deal with the effects of reducing the number of control volumes within the 
model and increasing the size of the trme increment. The number of control volumes is 
decreased from approximately fifty to only three. The size of the time step depends on 
the numerical method used and must be sized to provide stability in the formulation. An 
implicit method is used in the simple model and allows the time step to be increased from 
05 ms to 100 ms without loss of stability. Unfortunately, spatial resolution is decreased 
with an increase in the size of the control volumes and, if spatial resolution is important, 
the accuracy of the model will suffer. 

Specific control volume concerns include: 

a. Coolant Flow: Can coolant flow through the inlet plenum, particle bed, 
outlet plenum and the retention elements be simplified into a one dimensional flow com- 
pared to a two dimensional flow represented in the reference standard? 

b. Flow Distribution: Can simplifying assumptions be made to use a uniform 
flow distribution across the frits and particle bed? 

c. Sources of Pressure Loss: Can the sources of pressure loss in many small 
control volumes be combined in calculations for larger control volumes? 

d. Heat Transfer: Can the heat transfer prediction of the reference standard 


be well represented by using less control volumes? 
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Specific numerical concerns include whether the stability advantages of an implicit 
model are sufficient to counterbalance the longer computational times required. 

The following chapters describe the development and proposed application of the 
PBR and then describe the development of a less complex model to be used to analyze 


thermal-hydraulic transient response of a PBR fuel element. 


il 





CHAPTER 2 
BACKGROUND 


2.1 GENERAL 


The concept of operating a nuclear reactor system in space is not new. The require- 
ment to supply a reliable, sustained source of power to satellites and deep space probes 
was identified early in the US space program. Two programs were started to investigate 
the possible applications of nuclear reactors in space. The NERVA (Nuclear Engine for 
Rocket Vehicle Application) program was started in the early 1960’s and was primarily 
concerned with the development of nuclear technology for use as a means of propulsion 
in space. The SNAP (Space Nuclear Applications Program) program and Advanced Liq- 
uid Metal Cooled Reactor (ALMCR) program were developed in parallel with NERVA 
and were concerned with the development of power supplies for satellites and deep space 
probes (1960 through 1973). Although much of the SNAP program was restricted to the 
design and deployment of RTG (Radio-isotope Thermal Generator) power sources, an 
operating liquid metal reactor, the SNAP-10A, was launched in 1965 and subsequently 
tested at criticality in orbit. 

High power, pulsed, reactor concepts are being explored for orbital applications and 
general boost applications. With the renewed emphasis on space exploration, NTR (Nu- 
clear Thermal Rocket) technology is being reexamined for possible uses to minimize 
IMEO (Initial Mass in Earth Orbit) for long duration missions (L-2). Although the 
NERVA program has already established an NTR technology base, other advanced con- 


cepts are being studied as attractive alternatives to this 1960’s technology. 
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2.2 PARTICLE BED REACTORS 


2.2.1 Applications 


Placing any payload into LEO (Low Earth Orbit) requires extensive planning and 
support. This is usually measured by the cost of placing a unit weight into orbit (ie, $/kg) 
and places an emphasis on minimizing IMEO. Budget limitations drive the need for a 
compact, lightweight power system which can provide adequate energy for its applica- 
tion. Conventional terrestrial reactor designs are too bulky for use in space so advanced 
reactor concepts are required. One of these advanced concepts is the PBR which is 
projected to provide up to 50% greater specific umpulse, for a propulsion application, 
than the current NERVA design (L-2). 

PBR’s are being investigated by the Idaho National Engineering Laboratory (INEL) 
Multi-Megawatt Project Office, Brookhaven National Laboratory (BNL) and Sandia 
National Laboratory (SNL) (H-1, L-1, V-1). Open cycle systems (Figure 2.1) (M-1) as 
well as closed cycle systems (Figure 2.2) (G-1) are being examined. Open cycle systems 
are typical of a PEFNTR (pressure fed nuclear thermal rocket) (H-2) while closed cycle 
systems are typical of pulsed power applications being studied by SNL within the scope 
of the PIPE experiments being conducted in the Annular Core Research Reactor (ACRR) 
(V-1). PBR concepts show promise in that direct cooled particles within each fuel ele- 
ment provide a greater power density and a relatively large surface area for the heating of 


the working gas resulting in lighter and smaller systems (P-1,2). 
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Figure 2.1: Typical Open Cycle Nuclear Thermal Rocket System 
(Adapted from M-1) 
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Figure 2.2: Typical Closed Cycle Space Nuclear Power System 
(Adapted from G-1) 
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2.2.2 Particle Bed Reactor Design 


The previous investigation of PBR thermal-hydraulic characteristics was conducted 
by Tuddenham (T-1) and addressed the complications that arise during pulsed operations. 
Pulsed operations are applicable to open as well as closed cycle systems. A pulse asso- 
ciated with the open cycle system would be analogous to a timed burn in a chemical 
rocket system. 

The open cycle pulsed reactor usually consists of a coolant reservoir, a coolant 
pump, a preheating stage, the reactor assembly and the exhaust nozzle. These systems 
are discussed in detail in references B-2, L-2, P-1 and P-2. The coolant reservoir stores 
liquid hydrogen which is used as a coolant and eventually as the exhaust gas from the 
NTR nozzle. To maintain a relatively constant inlet pressure to the reactor assembly and 
provide sufficient mass flowrate through the core, a turbo-pump js placed in the system to 
force the coolant through the nozzle cooling jacket and then to the core. The hydrogen 
coolant is preheated becoming supercritical as it passes through the nozzle cooling jacket 
and innerpass cooling channels within the core. This provides efficient heat removal 
from the nozzle which is subjected to high temperature exhaust gases and ensures that the 
hydrogen coolant is in the gaseous state prior to entering the core. The hydrogen then 
passes through the fuel assemblies, is heated, and finally exhausted into space via the 


nozzle. 
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Figure 2.3: Packed Particle Bed Fuel Element 


(Adapted from L-3) 
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Figure 2.4: Packed Particle Bed Reactor Cross Section 


(Adapted from P-3) 
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As the hydrogen coolant passes through the fuel elements, it experiences several 
directional changes as shown in Figure 2.3 (L-3). The fuel elements are clustered, 
packed particle beds mounted in a moderator support assembly and surrounded by control 
drums and reflectors. A pressure vessel, most likely aluminum, encloses the entire 
assembly (P-2). The fuel elements may be clustered as shown in Figure 2.4 in a number 
necessary for the power required. 

Each fuel element assembly consists of fuel particles packed between the cold frit, 
the outer retention element, and the hot frit, the inner retention element. These retention 
elements are porous to allow the coolant to flow through the assembly but are fabricated 
to be robust enough to retain the fuel particles between them at elevated temperatures. 
The cold frit serves to distribute the coolant flow axially as well as to retain the fuel par- 
ticles. This allows better cooling of more power dense regions of the particle bed and 
prevent preferential flow through specific regions. The cold frit is made by sintering 
many small (2.5 um) stainless steel particles together. The hot frit, however, is made of a 
high temperature resistant material, usually rhenium, which has evenly spaced holes 
drilled to attain a desired porosity. 

The fuel particles are approximately 500 Um in diameter and consist of a central 
fuel kernel surrounded by two pyrographite layers and an outer layer of zirconium car- 
bide (Figure 2.5). Uranium carbide is used as the fuel and is enriched to approximately 
93% for compactness. These particles are packed directly in the annular region between 
the cold frit and the hot frit without being imbedded into a graphite matrix similar to high 
temperature gas reactors. Direct packing provides better heat transfer to the coolant and 
is expected to behave well during rapid power transients associated with pulsed power 


operations (P-3). 
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2.2.3 PIPE Experiment 


Pulsed Irradiation of a Packed Bed Element (PIPE) experiments have been con- 
ducted by Sandia National Laboratory (V-1). The PIPE experiments attempted to evalu- 
ate the performance of a packed particle bed fuel element in the areas of temperature 
characteristics, flow characteristics, fuel/coolant interaction, and power output. A fuel 
element test assembly (Figures 2.6 and 2.7) is then placed into the annular core research 
reactor (ACRR) for evaluation. 

Because previous modeling of a PBR fuel element concentrated on the specific 
geometry of the PIPE experiment, it is important to describe the experimental apparatus 
used for the PIPE experiments. The test assembly 1s a right cylinder 4.19 m long with a 
diameter of .35 m. Parahydrogen is introduced to the inlet plenum of the fuel element at 
a pressure of 2 MPa and a temperature of 300 K by a series of blowers which circulate 
the coolant through the test assembly. Flow enters the inlet plenum axially then flows 
radially inward through the cold frit, the particle bed and the hot frit. The hydrogen then 
exits the outlet plenum into a heat sink composed of stainless steel ball, through the flow 


meter and back to the blowers for another pass (V-1). 
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Figure 2.5: 500 Micron Fue! Particle 


(Adapted from V-1) 
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Figure 2.6: PIPE Experiment Fuel Element Test Assembly 
(Adapted from V-1) 
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Figure 2.7: PIPE Experiment Fuel Element Test Subassembly 
(Adapted from V-1) 
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Figure 2.8 shows detailed dimensions of the fuel element tested by Sandia National 


Laboratory. In addition to the measurements shown in Figure 2.6, other useful parame- 


ters are: 


¢Cold Frit Thickness 

eCold Frit Porosity 

eCold Frit Particle Diameter 
eCold Frit Material 

eHot Frit Thickness 

¢Hot Frit Porosity 


eHot Frit Material 
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Figure 2.8: Dimensions of PIPE Experiment Fuel Element 


(Adapted from T-1) 





2.2.4 Existing Models 


In his thesis (T-1), Tuddenham develops an analytical model to calculate the behav- 
ior of the PIPE experiment particle bed fuel element. His model, which will serve as the 
reference standard and which will henceforth be referred to as such, divides the fuel 
element into 50 control volumes. The reference standard is separated into 10 subdivi- 
sions in the radial direction and 5 subdivisions in the axial direction (Figure 2.9). Fur- 
ther, it uses a staggered grid arrangement (Figure 2.10) for discretization of the 
conservation laws. Table 2.1 specifies the geometry of the reference standard and shows 
how a varying cold frit thickness affects the volume of the inlet plenum and the particle 
bed. 

Although several other models exist for advanced space reactors (B-2, G-1, L-2), 
the reference standard addresses only the thermal-hydraulics of a PBR fuel element. It is 
more universal than one of the models used by Brookhaven National Laboratory (B-1), 
which assumes that power and coolant flow is matched in each control volume. The ref- 
erence standard separates these and provides detailed information pertaining to coolant 
temperature, pressure, and density for later use by a neutronic model. Also, the reference 
standard attempts to focus on actual operating conditions by using only instrumented sys- 
tem parameters, such as inlet and outlet temperatures and pressures, as sources of input. 

Tuddenham balances mass, momentum, and energy in each of his control volumes 
at each advance in time. A maximum timestep of 50 [Us is used by the reference standard. 
Although very detailed, the calculations tend to take a prolonged length of time. To 
examine a 6 second transient, for example, a single run would take approximately 48 
hours of computational time on an AT compatible desktop computer (80286). On the 


other hand, the detail provided in the reference model allows modeling of axial and radial 


26 





crossflows between the control volumes. Because the simple model combines many of 
the axial and radial contro] volumes into a single control volume, the particle bed for 
instance, specific information concerning crossflow is lost. If crossflow calculations are 


of primary importance, the reference standard should be used. 
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Figure 2.9: Control Volume Definition for the Reference Standard 
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Table 2.1: Reference Standard Control Volume Geometry 


(measurements in mm) 


Cold'Fut | | | Exit 
Axial Outer 
| Position 


Cold Fnit 
Thickness 


Inner PBED Hot Frit 
Thickness 


| Plenum 
Radius | Radius 
29.56 2.26 0.76 
29.39 1.95 0.76 
2950 1.79 0.76 
29,32 1.78 0.76 
29.38 1.91 0.76 


14.250 
14.250 
14.250 
14.250 
14.250 





(Volumes in Liters) 


ay ce 


Axial 
Position 


Height 
Volume 


53.00 
53.00 
53.00 
53.00 
53.00 


of Control | 


Volume 
of Inlet 
Plenum 


0.185 
0.187 
0.188 
0.188 
0.187 


Volume 
of 


Cold Fnt 


0.021 
0.018 
0.017 
0.017 
0.018 


hi Volume 


of 
PBED 


0.087 
0.088 
0.089 
0.089 
0.088 


Volume 
of 
Hot Fnt 


0.0037 
0.0037 
0.0037 
0.0037 
0.0037 


Volume 
of Exit 
Plenum 


0.0338 
0.0338 
0.0338 
0.0338 
0.0338 


265.00 | 0935 0.092 0.440 


TOTAL: 0.0185 0.169 


Void 
Fraction: 











Void 
Volume: 
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Figure 2.10: Reference Standard Staggered Grid Arrangement 
(Adapted from T-1) 
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CHAPTER 3 
THE SIMPLE MODEL 


3.1 MODEL BASIS 


A proposed particle bed reactor (PBR) which is to be placed in low earth orbit 
(LEO) will be subjected to multiple power transients. These transients may take the form 
of start-up and shutdown operations or rapid power transients during pulsed power opera- 
tions (B-1, G-1, L-1). Regardless of the form, any transient will require precise control 
of the reactor. As a part of the control system development, an accurate as well as a 
timely simulation must be accomplished via a digital model of the reactor. When cou- 
pled with an appropriate neutronic model, the reactor may be safely controlled by analyz- 
ing control module input parameters. Because fuel element neutronics are closely tied to 
the thermal-hydraulic behavior of the fuel element, it is important to develop a good 
thermal-hydraulic model. 

The basis of this research is to examine a method of analysis for the thermal- 
hydraulic response of a PBR fuel element in sufficient detail as to calculate reactivity 
feedback within the fuel for a subsequent neutronic model. The analysis should be, 
however, assiduously simple so as not to require excessive programming computations 
experienced with the reference standard model. Additionally, preliminary steps are to be 
taken to couple a single fuel element with other fuel elements within the core. 

The simple model will have some advantages over the reference standard, the most 
important of which is computational speed. Because the simple model is faster than the 


reference standard, many more fuel element variations may be investigated. This permits 


fot 





a wider scope of transient analyses and allows a departure from PIPE geometry. Further, 
the output data is presented quicker, granting the user a better understanding of cause- 
and-effect during the fuel element design phase. 

The following sections discuss the formulation of the simple model and the simpli- 
fications made to improve computational performance over the reference standard with 
respect to the amount of time necessary to perform a single analysis. Also in this chapter, 
the solid phase (fuel particles) and the gas phase (coolant) are described in detail. 

Finally, a method of solution is presented which combines the gas and solid phases into 


one model. 


3.2 SIMPLIFICATIONS 


One of the goals of the simplified model is to be able to model the thermal- 
hydraulic response of a fuel element in such a manner as to develop a real-time analysis. 
A real-time model is one which is able to assess a system’s response in a period of time 
less than or equal to the time it would take for an actual response. To do this, simplifica- 
tions must be made to the reference standard. 

Merely loading the reference standard code into a larger, faster computer (ie, Cray) 
would certainly reduce the amount of time necessary for a single run. However, the cost 
of analyzing the many variations would preclude the necessary translation. On the other 
hand, simplification of the mathematical method of the reference standard would produce 
a convenient program which could be executed on a typical desktop computer in a rea- 
sonable time frame. Further, a version of the model could be considered as part of a 


spacecraft power controller. 
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The two major simplifications incorporated into the simple model are the reduction 
of control volumes and the increasing of the time step. Reducing the number of control 
volumes provides a proportional decrease in the amount of computation required and 
increasing the time interval, reduces the number of times these calculations must be made 
during a specified simulation time interval. The overall effect is to reduce the net number 
of calculations. 

An example of the effect of the numerical simplification is readily seen when con- 
trasting the two models on a similar computer (80386). It takes approximately 8 hours to 
compute a 6 second simulated transient using the reference standard with a timestep of 
.O5 ms but only about 2 minutes for a 10 second simulation transient using the simple 
model with a timestep of 100 ms. This includes a reduction in the number of control vol- 
umes for the element from fifty in the reference standard to only three in the simple 
model. A computational advantage of approximately 240 to 1 1s realized with the simple 
model which allows analysis of 240 variations in the same time, on a comparable desktop 
computer, that it takes to do a single variation using the reference standard. 

The simple model uses only three control volumes compared to the fifty defined by 
the reference standard. The three control volumes are exhibited in Figures 3.1, 3.2 and 
3.3. The first control volume (Figure 3.1) includes the inlet plenum and the cold fit. 

The second control volume includes only the packed particle bed consisting of the fuel 
particles. Lastly, the third control volume includes the hot frit and the exit plenum. 
Combining the cold and hot frits with the inlet and exit plenums respectively allows a 
more precise model of the particle bed while defining conditions at the inlet and outlet of 


the fuel. 


33 





Stainless Steel 





Inlet Plenum 
Cold Frit 


a 
[| iSSSSSSSSESSSSSNESESIS 
® 





Control Volume fl: 
Inlet Plenum and Cold Frit 


Figure 3.1: Control Volume 1 
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Figure 3.2: Control Volume 2 
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Figure 3.3: Control Volume 3 
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3.3 SOLID PHASE 


3.3.1 General 


The solid phase of this model represents the fuel particles packed between the cold 
frit and the hot frit. For simplicity, the particle bed is modeled as a single control volume 
and, to simplify further, the particle bed is represented by a single, average fuel particle. 
Additionally, this average fuel particle 1s centered axially and radially in the particle bed. 
This approach is similar in nature to the one Tuddenham used for each control volume of 
the reference standard and assumes that all of the fuel particles behave in a analogous 
manner. 

Fuel particle dimensions and physical properties similar to the PIPE experiments 
are used in the modeling of the solid phase as shown in Figure 3.4. Once the physical 
properties of the fuel particle are defined, an overall heat transfer coefficient 1s developed 
from a weighted average of these properties and a surface fuel temperature may be deter- 
mined (M-2). 

Heat deposition rate also depends on the physical properties of the fuel particles 
used in the model. For a power transient, the particle bed power density is increased over 
a specific transient time to a predetermined level. The heat energy is not immediately 
transferred to the coolant however. The effects of heat deposition, heat storage, and heat 


removal must be taken into account. 


Sy 
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Figure 3.4: Fuel Particle Geometry 
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3.3.2 Energy Balance 


The fuel particles are the source of heat to the hydrogen coolant and can be modeled 
as a Separate control volume from the coolant. Since all of the fuel is contained in a 
single control volume, the following equation for the conservation of energy may be 


employed for the fuel: 


Soll 


= 
im 2) =q —hAWVcy(Ts — Tg) 


T = effective fuel particle temperature 


T’; = fuel particle surface temperature 

gq = heat source 

h = heat transfer coefficient 

Ay = particle surface area per unit volume 
Voy = size of the control volume 


Tg = bulk temperature of coolant 


Modeling of the solid phase for this model parallels the reference standard. Refer- 
ence M-2 uses a single node analysis for the particle bed. This allows the use of an 
assumed Steady state temperature distribution through the particle bed. In addition, 
material properties are assumed to remain constant. Reference M-2 simplifies the general 


heat balance equation in 3.1 by developing weighted average value: 


SF 


8/5) 


m,C, 


=q’’-U(T-Te) 
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m, = particle mass per unit surface area (kg/m’) 


C’, = average particle specific heat (J/kgK) 
q’’ = heat deposition per surface area (W/m’) 


U = effective over all heat transfer coefficient (W/m’K) 


It must be pointed out that T does not represent a temperature at a specific point in 


the fuel particle but is an average temperature for use in equation 3.2. The remaining 
values are determined in accordance with the procedures of M-2 which uses 1= 1 ... 4 to 


designated the shells of figure 3.4: 
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Although the properties associated with the fuel particle vary with temperature, the 
simple model, as well as the reference standard, hold these material properties constant. 
The exception is the heat transfer coefficient, 4, which is allowed to vary with tempera- 
ture and with coolant properties. A is shown in the following equation as a function of 
the Nusselt number and inversely with the sampling position in the particle bed. This 
relationship will be discussed in greater detail in the next section. 


3.4 


Keeeiard 


x 


h Nu 


coolant ~ P 


where: x = position in the particle bed (in this case x will be half of the particle bed 
thickness) 


Conductive heat transfer between adjacent fuel particles and control volume bound- 
aries as well as radiative heat transfer, which can occur at elevated temperatures, are 
neglected in the simple model. These forms of heat transfer are addressed and discussed 
at length in reference T-1 which concludes that their effects are minor compared to the 


convective heat transfer. 


3.3.3 Discretized Form 


In order to analyze the thermal-hydraulic response of the PBR fuel element, numer- 
ical approximations must be made to the ordinary differential equations and the convec- 
tive equations used in the model. The solid phase heat balance equation addressed in 
equation 3.2 may be solved by using a discretized form of the equation. 

As with the reference standard, each side of equation 3.2 is multiplied by the con- 
trol volume size, Vy, and the particle surface area per unit volume, Ay, to produce total 


power terms (Watts). This is shown in the following expression: 


4] 





= dT = = BES 
m,C pVcv4Ay— =q’ VaAy —- UV pyA(T -— Ta) 


where Ay can be expressed as a function of the void fraction and the fuel particle 


diameter, D>: 


_ 6(1-€) 3.6 


A 
V Dp 


To simplify this expression, incorporate equations 3.7, 3.8, and 3.9: 


Let A =UV Ay 3.7 
Let B=m,C ,VcvAy 3.8 
Let O = q ’V cyAy 3.9 


These substitutions in equation 3.5 result in: 


pT _ 3.10 
Bo = OAT =i) 


A simple implicit scheme is used to place equation 3.10 into a discretized form. 
This helps to provide stability at larger timesteps. In addition to using T at time n+1, Q 


n+1 


is also advanced ton+1. Q”” represents the driving element in this relationship. If Q 
were to remain unchanged, the result would be a null transient. Further, because the A 
(equation 3.7) is a function of properties which change with temperature, it 1s evaluated at 


time n for the discretized form of 3.10 which can be expressed as: 
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When this is solved for T””', a simple expression results: 


SeleZ 


—n+]1 


oi mC ae 
P 1464 
Equation 3.14 is incorporated into the simple model to advance the average temper- 
ature of the fuel. The interphase heat energy transferred from the solid phase to the gas 
phase can now be calculated from the temperature difference between the average fuel 


temperature and the bulk temperature of the coolant. The interphase heat transfer, Q,, is 


expressed as: 


of =04,Vo(7""—11) 3.13 
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3.4 GAS PHASE 


3.4.1 General 


Although less complicated in nature than the reference standard, many of the fea- 
tures incorporated in modeling the gas phase in the reference standard are included in the 
simple model. The gas phase is, however, more complex than the solid phase in that it 
includes both hydraulic and thermal characteristics. Complications arise, however, in 
defining the thermal-hydraulic balances within each control volume. 

The general structure of the gas phase equations balances energy, mass, and 
momentum for each control volume. As enthalpy and mass flowrate are solved for each 
control volume, these values are advanced in time and forwarded to the adjoining control 
volume. This allows the model to generate a time response for the fuel element during a 
simulated transient. 

With the exception of the total number of control volumes, the simple model 1s 
based on the same assumptions and significant features which affect the element response 
addressed in the reference standard. As in the reference standard, the simple model is 
cylindrical and includes axial and radial flows. Cross axial flow in the particle bed is not 
a variable for the simple model since the flow vectors are contained entirely within the 
second control volume. 


In the previous section, it was mentioned that the heat transfer coefficient, h, for the 


hydrogen coolant varies with temperature and velocity. Reference E-1 provides an 
expression for the heat transfer coefficient in terms of the Nusselt number, the heat con- 


ductivity of hydrogen and the sample position within the particle bed: 
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where: 


Nuk 3.14 





Nu, = Nusselt number; 


k = coolant heat conductivity; 


xX = .5 (particle bed thickness) 


Reference E-1 also provides an expression for the Nusselt number: 


where: 





Nu, =0.8Re’Pr> QMS 
vd 
Re= lial 
UL 
1-8 
> 
6(1 —€) 
aa 


P 


D,, = particle diameter; 


d,, = effective particle diameter (E-1) 


The rest of this section discusses the energy, mass, and momentum balance tech- 


niques and equations used in the simple model. Once the basic formulations are 


identified, a discretized form of the relationships will be stated and terms specific to each 


control volume will be assembled for final analysis. 
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3.4.2 Energy and Mass Balance 


The primary energy transfer mechanism within the particle bed is convection heat 
transfer and it is assumed that all of the heat deposited in the fuel particles will be trans- 
ferred to the hydrogen coolant. In order to accurately calculate the response of the entire 
fuel element to changes in input heat energy, an energy balance and mass balance must 
be performed on each control volume. Since there are fewer control volumes in the sim- 
ple model as compared to the reference standard, the calculations involved will be fewer 
as well. 

The interphase heat transfer provided by equation 3.13 represents the heat energy 
being convected to the hydrogen coolant. This serves to couple the solid phase with the 
gas phase of the model. The rate at which heat is transferred from the solid phase to the 
gas phase depends on the rate of increase of fuel heat deposition and coolant properties. 

For a control volume with fixed volume, the time rate of change of enthalpy is the 
sum of the heat energy source, the time rate of change of control volume pressure and 
volume, and the sum of enthalpy flux terms. This is similar to equation 3.22 in reference 


T-1: 





d(MH) dp i! eG 
— Vo — + W.H. 
( at 2:4 Y dt 2 % 


— 


where: M = mass within the control volume; 
Q, = heat transfer from solid to coolant for control volume; 
P = average pressure within the control volume; 
H = specific enthalpy; 


W = mass flowrate. 
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For the simple model fewer control volumes limit the number of enthalpy flux 


terms and, as a result, equation 3.16 may be written as: 


d(MH dP SF 
=O) + Vey + W, H. Ano ieee 


in in 





The general mass balance equation may be expressed as: 


dM 3.18 
= W., a ee 
dt 


Equations 3.17 and 3.18 may be combined to produce a more direct equation which 
will allow us to advance enthalpy. This is accomplished by expanding the left hand side 
of equation 3.17, multiplying both sides of equation 3.18 by enthalpy, H, and subtracting 
equation 3.18 from 3.17. Also assuming that H,,, 1s representative of the enthalpy in the 
control volume (a donor cell method) and mass flowrate, W, is the same throughout the 


control volume, a general relationship 1s produced: 


A ou, 7 aP ENG ees 3.19 
dt mn” CV dt Q, ( in fe) 





In the discretized form of equation 3.19, and adopting an implicit formulation, 
values at the current timestep may be used to compute the enthalpy for the subsequent 
timestep. An implicit numerical formulation is also used with the enthalpy to provide 


numerical stability at large tumesteps. This results in: 
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Values of Q;'*' and H;,*' are determined from the previous control volume (or inlet 


boundary) and may be evaluated at n +1. For example, H/,*' for the first control volume 
represents the enthalpy of the coolant being supplied to the control volume. The source 
code allows a change in the inlet temperature during a transient which would subse- 


quently affect the enthalpy of the coolant leaving the control volume. 


3.4.3 Momentum Balance 


In addition to a balance of mass and energy, a balance of momentum must be dove- 
tailed into the response of the fuel element. For the simple model, the framework of the 
MINET (Momentum Integral Network) code (V-2) is used to determine the effects of 
momentum on the flow of coolant. The MINET code develops equivalent forces which 
act to resist the flow of the coolant through each control volume. These resistive forces 
can be divided into the significant mechanisms which contribute to the head loss through 
a fuel element. Unlike the MINET code, which accounts for spatial changes in mass 
flowrate, the simple model assumes that mass flowrate is the same everywhere (at any 
time t). 

The basic momentum integral equation (V-2) used by the simple model is 


dw Salk 


lov 


N 
= lap = ete ~ C5) 
i=] 
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where: /¢,=L/A (control volume inertia) 


L = average length of travel for each control volume 
A = cross sectional area perpendicular to flow 


F’., = equivalent resistive pressure drop 


This equation is equivalent to a constitutive relationship which relates a change in 


pressure to the flowrate through a control volume (C-1): 


AP = RrotalW 3.22 


where: Rr,,,; = total equivalent resistance; 


AP = pressure drop. 


Equation 3.21 may now be expressed as a function of the total pressure drop across 
the control volume and the sum of flow resistances: 


aw 51045: 


oe = AP cy — D(Ry)W? 
This is a convenient expression because all of the terms may be determined from 
existing relationships. These include relationships for an equivalent pressure drop due to 
friction, spatial acceleration, expansion or contraction, manifold mixing, and packed 
particle beds. Although the MINET code includes gravity terms, they are considered 
negligible. Expressions which differ from the reference standard will be examined more 


closely. 
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The equivalent pressure drop due to friction is calculated for the inlet and the outlet 


plenums only and takes the form of equation 3.24: 


F=f pe al 3.24 
if jie 2pA? 


where: F,= equivalent pressure drop due to friction; 








L, = plenum half-length; 
fp 0leeRe '' ; 








R pvD., 8) 5 
Cc — 
UL 
4A 
on 


Since frictional pressure drop is directly related to the distance that the coolant must 
travel through the control volume, an average path length is defined for the inlet and 
outlet plenum which is half of the overall element length. The cross-sectional area used 
in this calculation is the area in the plane perpendicular to the axial axis of the fuel 
element. This area is used because the path length in the axial direction (for the inlet and 
outlet plenums) is so much greater than the radial path length. Density is the average of 
the density calculated at the inlet and outlet of each plenum. 

Equation 3.25 relates a correlation line which was fit to the log-log Moody plot of 
friction factors as a function of Reynold’s number, Re, and roughness in the channel 


caused by the frits (T-1). This improves on the relationships developed by McAdams and 
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Blasius to correlate friction factors for smooth pipe and is the relationship used by the 
reference standard. Resistance due to friction is determined in the inlet plenum and outlet 
plenum only since the Ergun (E-2) relation includes this effect for the particle bed. 
Spatial acceleration in a duct which allows cross sectional area to change while 
holding coolant density constant will produce a change in pressure across a control vol- 


ume. These conditions are experienced in the inlet and outlet plenums and is: 


z i hie 3.26 
3S ARE AG GS 


where: F’, = equivalent pressure drop due to acceleration. 


The cross-sectional areas in this case are the areas associated with the inlet and 
outlet of the plenum control volumes. For example, for the outlet plenum, A, is the area 
at the exit of the fuel element and A, is the cross-sectional area at the inlet of the hot frit. 
A, for in the inlet plenum is the cross-sectional area at the outlet of the cold frit and A, is 
the area at the inlet to the fuel element. As in the calculation for the frictional pressure 
drop, the density used in equation 3.26 is the average of the density at the inlet and outlet 
of the control volume in question. 

Equation 3.26 can not be used for the control volume containing the fuel particles 
because the coolant experiences a change in duct area and density. A general form of the 
change in pressure caused by spatial acceleration in a duct which allows cross sectional 
area and density to change is addressed in Appendix A. This case describes the changes 
to the hydrogen coolant as it travels radially through the particle bed. The coolant wil 
experience a change in cross sectional flow area as well as a change in density due to the 


transfer of heat. This takes the form of: 


all 











A, +A, A, +A, 3.27 
es a i 2A,A>, Wy, _ Da Wy, 


where: (F,) pac) = Spatial acceleration losses in the particle bed; 


A, = cross-sectional flow area at the inlet of the particle bed; 
A, = cross-sectional flow area at the outlet of the particle bed; 
v, = coolant velocity at the inlet of the particle bed; 

v, = coolant velocity at the outlet of the particle bed. 


Another source of resistance within the fuel element is the sudden expansions and 


contractions experienced by the coolant as it passes through the cold frit and hot frit. 
These expansion/contractions would occur when the coolant enters and exits each frit 
because of the area change on both sides of the frits. For example, the coolant would 
experience a first contraction as it entered the cold frit and a second contraction as It 
leaves the cold frit and enters the particle bed. The coolant would experience a first 
expansion as it leaves the particle bed and passes into the hot frit and a second expansion 
as it enters the outlet plenum. The general form of this equivalent pressure drop is 


adopted from T-2 and is given as: 


Ww? 3.28 
1G 5 
2PA Snall 


i= 





where: F, = equivalent loss due to a sudden expansion or contraction; 


K,=K, for expansions; 
K,=(1-B)’; 


K,=K,, for contractions; 
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A mall = the smaller of the two areas involved. 

The value of density used in equation 3.28 is the density at the point of expansion or 
contraction. In the case of the cold frit, the density calculated at the entrance of the 
particle bed is used to evaluate the equivalent pressure drop due to sudden contractions. 
The density calculated at the exit of the particle bed, on the other hand, is used for the 
expansions observed at the hot frit. 

The final relationship to be used to describe the hydraulics of the fuel element is for 
the mixing effects experienced in the inlet and outlet plenums. This is known as man- 
ifold mixing. The general relationship was derived by Bajura (B-4, 5, 6) and improved 
by Datta and Majumdar (D-1, M-3), and is used in the reference standard. When put into 


a form compatible with equation 3.22, the manifold effect 1s expressed as: 


Fry = OXF y —adiat + Fy — axial) 6 Ga, 
Ww? 
he aia 
‘i Ae ial 
Ww? 
Fy radial om La 
PA padial 


where: Fy, = equivalent pressure drop due to the manifold effect 


© = .95 for the inlet plenum and 1.1 for the exit plenum 
The value of density used in this case is, again, the average of the inlet and outlet 


densities for the control volume being analyzed. In addition, it must be noted that 


although the reference standard uses a value of 2.66 for exit plenum © in the calculations, 
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T-1 states that the values of © may be adjusted to take into account differences in mixing. 
Because of the difference in how the simple model treats crossflow within a control 
volume, © for the outlet plenum is reduced from 2.66 to 1.1. This allows a better 
approximation of the results obtained by the reference standard. Although this is minor 
contributor to the total flow resistance, it serves to illustrate how the model may be fine 
tuned and results in a better agreement for pressure drop-flow relations. 

To complete the set of expression used in the simple model, the Ergun relation is 
used to determine the equivalent pressure drop across the cold frit and the particle bed. 
Since the cold frit is made of small particles, modeling it as a particle bed is a good 
approximation. The Ergun relation treats the flow through a packed particle bed as flow 
through tubes with irregular cross sections vice flow around many objects and is the gen- 
erally accepted approach to modeling particle beds (B-6, E-2). The general form of the 


Ergun relation is: 


s 150uy,(1 2} - 1.75pv,1¥,]( 3 3.30 
Phed — +b rae b 
Die D,€ 


where: Fp,,, = equivalent pressure drop due particle bed effects; 


L, = thickness of the particle bed; 
D,, = particle diameter; 

V, = superficial velocity; 

Lt = viscosity; 


€ = void fraction or porosity; 
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Using the conservation of mass relation, this expression may be put into a form 


which is compatible with equation 3.27 and will appear as: 


1500ntl—cy 7sl—< 3.31 
Fn Du "sf ( a 
DE DAW DEE DA, 





Once more, the density used in this relationship an average of the density at the 
particle bed inlet and the density at the exit of the particle bed. 

Having defined all of the significant terms used to describe the pressure loss 
through the fuel element, the next step is to discretize equation 3.21. An implicit scheme 
may be used if W’ is separated into a combination of the mass flowrate evaluated at n and 


n+1 


n+1 to give:W"W""’. The discretized form of equation 3.21 is then: 


(W" +1 ~Ww”) 3.32 


J n n n n 
St = 7 (OP cy — Rona W W ") 
cv 
Solving for W”*’: 


n+1 


W"Icy + StAP2, 3.33 
Tey + 8tRe iia W" 
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This section has described the process by which the simple model is able to advance 
enthalpy and mass flowrate. The following section will describe how this is coupled with 
the advance of fuel temperature and interphase heat transfer to produce a full portrayal of 


the thermal-hydraulic response of a typical fuel element. 
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3.5 METHOD OF SOLUTION 


A simple flow chart of the transient solution is shown in Figure 3.5. As can be seen 
in the schematic, initial conditions are established by the user and an initial, steady state 
condition is calculated. This provides initial state variables to the transient portion of the 
source code. Next, enthalpy and mass flowrate 1s advanced for each control volume. In 
addition, in control volume 2, calculations include the advance of fuel temperature. This 
in turn provides the calculated value of the interphase heat transfer. 

Of note, the control variables (heat deposition, inlet pressure, outlet pressure, and 
inlet temperature) are also advanced in accordance with the transient parameters selected 
at the initiation of the simulation. 

Coolant properties, such as density, viscosity and heat transfer coefficient, are recal- 
culated for each control volume based on new temperatures and pressures corresponding 
to the value of enthalpy and mass flowrate. Finally, fuel element information is sent to 
the output file and time is advanced by the chosen timestep for the next set of calcula- 


tions. 
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Figure 3.5: Flow Chart for Transient Solution 
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CHAPTER 4 
VALIDATION AND APPLICATIONS 


4.1 GENERAL 


Although originally intended as an interim step toward the development of a real- 
time controller for a particle bed reactor, the simple model also lends itself serendipi- 
tously as a design tool. Because the time necessary for a single analysis is significantly 
shorter than that required by the reference standard, many more variations may be 
examined and compared. Prior to using the source code to examine these different varia- 
tions, the source code must first be verified as a valid model. This may be done by com- 
paring the results from the simple model with the reference standard. 

An additional feature of the simple model is that it is capable of analyzing transients 
in which control variables change (inlet pressure, outlet pressure, inlet temperature, and 
heat deposition rate). That is, inlet pressure, for example, may be permitted to increase or 


decrease during the transient. 


4.2 VALIDATION 


Validation of the simple model is accomplished by comparing results of steady state 
and transient calculations provided by each model. Also, as a verification that the numer- 
ical scheme is valid, null transients are analyzed at three different power density levels. 
Overall, the simple model agrees very closely with the reference standard but lags 
slightly during a baseline 1 second transient from 0 GW/m’ to 2 GW/m’. Null transients 


are performed at 0, 1, and 2 GW/m’. 
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A more in depth discussion of the validation procedure is attached as Appendix C 
and includes graphical comparisons of the simple model and the reference standard. 

After completing a successful validation, the simple model may be used to examine 
variations other than the baseline (PIPE) configuration used by the reference standard and 


the PIPE experiment. Examples of which are supplied in the following sections. 


4.3 DESIGN VARIATIONS IN FUEL ELEMENT GEOMETRY 


The source code for analyzing the response of the packed particle bed fuel element 
allows changes in the dimensions of the fuel element. Geometry variables which may be 
changed include outer element radius, outer cold frit radius, cold frit thickness, particle 
bed thickness, hot frit thickness, and length. Other parameters which may be changed 
which influence the hydrodynamic characteristic of the fuel element include particle 
diameter, frit porosity, manifold mixing coefficients, and fuel void fraction. 

As a demonstration, mass flowrates for a range of fuel element power densities 
were determined for a number of different fuel element lengths and plotted against the 
thermal power that they would deliver. Figure 4.1 shows the results of these calculations 
for lengths of .265 m (baseline), .500 m, .750 m, and 1.000 m. In each case, inlet pres- 
sure was maintained at 2000 kPa; outlet pressure was maintained at 1915 kPa; and inlet 
temperature was maintained at 300 K. Also, each curve has eleven points of data which 
represents calculations based on different power densities ranging from 0 GW/m‘’ (ther- 
mal power equal zero) to 1 GW/m’ (the end point of each curve) at .1 GW/m’ intervals. 

For an open-cycle nuclear thermal rocket application, Figure 4.1 could aide the 
designer chose a range of power operations given specific mass flowrate restrictions. Or, 


on the other hand, the designer could calculate the amount of hydrogen fuel/coolant 
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required given power level and fuel element size restraints. Figure 4.1 could be paired 


with a plot of exit temperature verses thermal power to ensure component temperature 


limits are not exceeded. 
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Figure 4.1: Mass Flowrate vs Power for Various Element Lengths 
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4.4 VARIATIONS IN THE MANNER OF TRANSIENT CONTROL 


The code has the provisions to fix the initial and final values of inlet and outlet 
pressure and inlet temperature at a single value. They may also be changed linearly dur- 
ing a portion of a transient of specified duration. The time during which one of these 
variables may change will typically start when the power transient begins but 1s not 
coupled to the same transient duration. For example, power density may be increased 
from 0 GW/m° to 1 GW/m in 1 second and inlet pressure decreased 90 kPa over a period 
of 5 seconds. To conserve fuel, a nuclear thermal rocket may be operated in such a fash- 
ion, that is, the reactor 1s started up prior to admitting the hydrogen coolant through the 
fuel element. Figures 4.2 and 4.3 illustrate this operating mode for a 1 m long fuel 
element and show how mass flowrate and delivered thermal power (the product of mass 
flowrate and enthalpy change) respond to such a transient. 

Of note, Figure 4.2 indicates that thermal power will overshoot the required neutron 
power due to the extended decrease in the outlet pressure. While an increase in neutron 
power tends to increase the outlet temperature, a decrease in outlet pressure tends to 
increase mass flowrate. Thermal power continues to increase until outlet pressure 


reaches the required operating value of 1900 kPa. 
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Thermal Power Response to Combined 
1s Power/15s Outlet Pressure Transient 
(L=1m; P-Inlet=2000 kPa; T-Inlet=300 K) 
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Figure 4.2: Thermal Power Response to Combined Transients 
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Mass Flowrate Response to Combination 
1s Power/15s Outlet Pressure Transient 
(L=1m; P-Inlet=2000 kPa; T-Inlet=300 K) 


Flowrate (g/s) 


60.0 
50.0 
40.0 
30.0 
20.0 
0.0 4.0 8.0 12.0 16.0 20.0 24.0 
Time (s) 


Transient Starts at t=ls 


Figure 4.3: Mass Flowrate Response to Combined Transients 


64 





4.5 CORE REPRESENTATION BY MULTIPLE FUEL ELEMENTS 


As a further demonstration of how the simple model may be used as a design tool, a 
cluster of 19 fuel elements are arranged into three rings around a center fuel element as 
presented in Figure 4.4. The center element is designated as O and is surrounded by the 
A-ring, B-ring, and C-ring. Using a zero order Bessel function of the first kind to 
approximate the relative power densities of each ring, the simple model was used to cal- 
culate the time behavior of each ring of the reactor for mass flowrate, thermal power, and 
exit temperature. Length was set to 1 m, inlet pressure to 2000 kPa, outlet pressure to 
1900 kPa, and inlet temperature to 300 K. 

Figure 4.5 shows the mass flowrate response to a 4 second power transient for a 
single fuel element in each ring. The heat deposition rate is increased from zero to 1 
GW/m’ in element O with the other rings being increased proportionally. As expected, 
the fuel elements with the lower power densities develop less flow resistance and, as a 
result, mass flowrate 1s greater through these elements at the final power level. Mass 
flowrates specific to each element may be combined to describe the overall mass flowrate 
response to the transient as shown in Figure 4.6 such that: 


19 4.1 


Figure 4.7 depicts total thermal power response to the 4 second increase in neutron 
power. The total thermal power is found by summing the product of mass flowrate and 
change in enthalpy for each element: 


19 4.2 
> WAH, 
p=] 


i= 
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Finally, Figure 4.8 describes exit temperatures for each ring of elements. The 
average exit temperature represents the temperature of the coolant which is delivered to 
downstream power plant components and is determined by weighting each exit tempera- 
ture with the corresponding mass flowrate. Note that the center ring will be operated well 
above the average outlet temperature which indicates that a limit associated with fuel 
element material may be reached before any other component (ie, turbine blades). 

As a design consideration factor, outlet temperatures may be made more uniform by 
changing the orificing of the inlet plenum or the thickness of the cold frit for the A, B, 
and C-rings. This would serve to decrease thermal gradients across regions of the core as 


well as improving allowable thermal power. 
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Figure 4.4: Clustered Fuel Element Armangement for 19 Element PBR 
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Figure 4.5: Mass Flowrate Response for 19 Element PBR 
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Figure 4.6: Total Mass Flowrate Response for 19 Element PBR 
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Thermal Output Power Response to 4s 
Up-Power Transient for 19 Element 
Particle Bed Reactor 
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Figure 4.7: Total Power Response for 19 Element PBR 
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Outlet Temperature Response to 4s 
Up-Power Transient for 19 Element 
Particle Bed Reactor 
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Figure 4.8: Exit Temperature Response for 19 Element PBR 
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4.6 REAL-TIME CONTROLLER 


The simple model takes a step toward development of a real-time model for incor- 
poration in control of a space-based reactor by significantly decreasing the time necessary 
for computation. Although some information is lost when compared to the data delivered 
by the reference standard, the overall framework which describes fuel element perform- 
ance is similar. 

As a demonstration of real-time control, the simple model was modified in a man- 
ner which suspended the process of writing data to a data file. The calculation time for a 
transient was less than the transient simulated (for a single element). A 10 second 
transient, for example, required approximately 3 seconds of computation time for a time- 
step of 100 ms and approximately 7.5 seconds for a timestep of 10 ms. The calculations 
were performed on a desktop computer equipped with a 30386 microprocessor (16 MHz). 

Although there are reasons that these results are inappropriate for a final control 
application, they do indicate the plausibility of using a version of this thermal-hydraulic 


model. The reasons for the continued adjustment of these results are: 


¢ Improvements in computer speed are expected; 


¢ An optimization of calculation techniques has not been performed; and 


¢ Multiple elements and other computations must be done in parallel. 
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CHAPTER 5 
CONCLUSIONS AND RECOMMENDATIONS 


5.1 CONCLUSIONS 


The previous sections describe an investigation to construct a simple model which 
could effectively analyze the thermal-hydraulic response of a packed particle bed reactor 
fuel element during steady state and transient conditions. The fuel element was first 
divided into a solid phase and a gas phase which are modeled separately then coupled for 
the final solution. 

The solid phase 1s modeled in a single control volume and a lumped parameter 
approach is used. This approach is essentially the same as the approach used in the refer- 
ence standard and generates a fuel temperature which is representative of an average fuel 
temperature instead of a local fuel centerline temperature for instance. Interphase heat 
transfer is calculated from the product of an overall heat transfer coefficient and the dif- 
ference between the average fuel temperature and the bulk temperature of the coolant. 
This interphase heat transfer rate provides the heat input to the coolant as it passes 
through the fuel particles. 

The gas phase is modeled in three control volumes: the inlet plenum combined with 
the cold frit; the packed particle bed; and the outlet plenum combined with the hot frit. 
This particular arrangement allows the fuel region to be modeled separately from other 
components of the fuel element. 

Conservation equations for the conservation of energy, mass, and momentum (in 


the form of the momentum integral equation) are used to ensure a balance when analyz- 
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ing the fuel element during a transient. Steady state values and transient values (for a1 s 
baseline transient from 0 to 2 GW/m’) are analyzed and then compared to the reference 
standard. 

The calculations resulting from the simple model compare favorably with the refer- 
ence standard and, because of the simplifications, are able to analyze the same transient 
significantly faster as well as compute steady state values directly. After eliminating 
writing to a data file, the simple model is able to achieve a faster-than-real-time thermal- 
hydraulic analysis of a transient for timesteps larger than 10 ms and an 80386 micropro- 
cessor. 

The simple model may also be used as a design tool. It may be used on a desktop 
computer and many different calculations may be accomplished during a reasonable 
timeframe. Further, the simple model has the capability to accept variations in geometry, 
fuel composition, and initial/final conditions. 


5.2 RECOMMENDATIONS FOR FURTHER STUDY 


Development of the simple model is another step in the creation of a real-time, 
autonomous controller for a packed particle bed reactor. Although this model represents 
the thermal-hydraulic portion of the controller, it sets the stage for further investigation. 

Some areas which would merit further investigation include: 

¢ Calculations obtained by the reference model and the simple model should be 
compared to experimental results. At the time this investigation concluded, such exper- 
imental results were not available. 

¢ A neutonics module needs to be developed which uses the data generated by the 


thermal-hydraulics module. A real-time controller may then be constructed and tested. 


74 





¢ A model which includes additional power generating system components should 
be developed. Components, such as a high temperature turbine, turbo-pump, and asso- 
ciated valves and piping could be coupled to the representation of the core provided by 
the simple model. 

¢ Study a means to adjust the simple model to correct the observed sluggishness 
(Section C.3) 

¢ Continue to compare the results obtained by the simple model with future models 


which may construed as a reference standard. 
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APPENDIX A: 
DERIVATION OF THE GENERAL CASE FOR SPATIAL ACCEL- 
ERATION IN A DUCT WITH 
NON-CONSTANT AREA AND NON-CONSTANT DENSITY 


A.1 BACKGROUND 

Spatial acceleration term in a momentum balance equation are usually presented as 
the pressure change in a duct caused by a change in density or by a change in cross sec- 
tional area. In the former case, the pressure change in a duct caused by a change in den- 


sity while maintaining cross sectional area constant is expressed as: 


P, —P,= PV — PV; A.1 


In the latter case, a change in pressure associated with spatial acceleration in a duct 


with constant density but with a change in duct area, is expressed as: 


leh \ ye A.2 
where W = Mass Flowrate 
In the case of the packed particle bed reactor fuel element, the working gas is 
subjected to a combination of cross sectional area change and density change as it travels 
in the radial direction. Although there are integration methods available which could 
accurately describe this fluid behavior, it would necessitate dividing the particle bed into 


many control volumes. For the simple model, it 1s desirable to maintain the particle bed 
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within a single control volume. Because neither equations A.1 nor A.2 accurately 
describe the thermal-hydraulics associated with variable density gas through the particle 
bed, a general case derivation must be formulated for the simplified model. 
A.2 THE GENERAL CASE 

Figure A.1 shows the general case in which the cross sectional area in the duct and 
the density of the gas is allowed to change as the gas travels through the duct. The objec- 
tive of this derivation 1s to formulate an expression which will satisfy the limiting cases 
described in equations A.1 and A.2 above. To do this we must first define an average 


pressure P as acombination of P, and P,: 


P =aP,+(1—a)P, A.3 


Where a and (1 — a) represent fractional constants between 0 and 1. 


W in— Wout 


Figure A.1: Duct with changing area and changing density. 
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To begin the derivation, a momentum balance analysis must be performed on the 
control volume shown in Figure A.1. P becomes important at this point because it is 
allowed to act at some point within the duct on an area which is defined as the difference 
between A, and A, when projected onto the vertical plane. This then results in the 


momentum balance equation: 
P,A,— PA, + P(A, —A,) = PovzAy — piviA, A.4 
Substituting the definition of P into equation A.4 and collecting terms gives: 


(P,—P,){A,+a(A,—-A))} = 92V3A, z VIA, oD 


By inspection, equation A.5 reduces to equation A.1 when A, and A, are set equal 
to one another. When the duct inlet area and outlet areas remain the same, the term 
a(A, —A,) is eliminated resulting in the limiting case shown in equation A.1. For the 
other case in which the solution is known, the fractional constant @ 1s solved while 
density remains constant. 

For the case in which density 1s held constant, the Bernolli relationship 1s used to 


find an expression for the pressure change on the left hand side of equation A.5: 


A.6 
P,-P,=5 03 -¥) 


and, because mass is conserved in a steady state condition, constant cis used to demote a 


ratio of areas and hence a ratio of velocities: 
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| A.7 


Substituting the expression for P,—/P, and letting v,=cv, results in the following 


expression: 


Joes A.8 
ale —1){A, +a(A, -A,)} =A,(c — 1) 


The fractional constant a can now be solved by dividing both sides by A,(c — 1): 


(c+1){c +a(1-—c)}=2c 
c(c+1)+a(c+1)U-—c)=2c AY 
a(c+1)(c -1l)=c(c-1) 


and upon completion of the algebra, an expression is derived for the fractional constant 


a: 











aS 





Equations A.10 and A.11 are then substituted into equation A.3 for an expression 


= A, A, A.12 
P= ,+ 2 


Finally, substituting this new expression for P into equation A.4 and solving for the 


for P: 








pressure change provides an expression for the general case: 


> —_p _( Ait4e (Arts | A.13 
AAD 2A 
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APPENDIX B: 
SUMMARY OF EQUATION FOR HYDRAULIC ANALYSIS 
OF A PBR FUEL ELEMENT 


1. Constitutive Relationship: 
AP, = R,;W" 
AP, = Total pressure drop across element 


R, = A term representing the total resistance to fluid flow 


W = Mass flowrate through the element 


2. Momentum Integral Equation: 
MINET CODE (Reference V-2) 


dw 
I —=P.-P,+F 
m dt i oO Ji 


ressEquivalent 


ia 


PressEquivalent 


=LUF,+F,+F,+F,+F,+F;) 


a 


PressEquivalent 


es 
=L(R, +R, +R, +R, +R, +RW 


Rp=X(R, +R, +R, +R, +R, +R,) 
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3. Pressure Relationships for Inlet and Outlet Plenums: 


L WwW? 
Dig J2PA 


f,=0.138Re~ 





L, = The half-length of the plenum 


1 1 \w 
a ie aaa 
A; A; }2p 


Ww? 
et ar 
PA DENS 
K, = K for expansions K,=(1-BYy 
K, = K_forcontractions l 
a K,==(1-By 
Z 
_ smallerarea 
~ largerarea 


Fy = Oy rer seo) eee, 


2 
W 
M -—axial ~ 9) 
PA Zia! 
and 
WwW? 
18 M ~radial — A? 
p radial 


© = .95 for Inlet Plenum; 1.10 for Outlet Plenum 
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4. Pressure Relationships for Packed Particle Bed (also Cold Frit): 


(Ergun Relationship) 


150uv,(1 —)- 1.75pv,|v|(l—-e 
Lv ( Th eu pv, iv I : 


F Phed — Lj D BE iB, 


OY: 


150u(1 —¢) 1.75(1-¢) \|_, 
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APPENDIX C: 
MODEL VALIDATION 


C.l STEADY STATE 


C.1.1 General 

In order to provide one set of validations for the simple model, it may be compared 
to the reference standard. Steady state and transient results for both models are compared 
graphically in the following subsections. Values for reference standard parameters were 
drawn directly from T-1 and values for the simple model were obtained from the steady 
State and transient source codes attached as Appendices C and D. Also, each model used 
the Sandia National Laboratory PIPE experiment fuel element as the basis for its overall 


configuration. 


C.1.2 Mass Flowrate 

At a Steady state condition, coolant mass flowrate will vary as a function of the par- 
ticle bed power density. That is, the mass flowrate will adjust itself according to the 
resistance to flow produced by the element as a whole which will be a function of 
temperature and pressure. Figure C.1 shows how mass flowrate varies as a function of 
power density. Keep in mind that the average temperature of the fuel and the exit tem- 
perature increase with an increase in power density. Also shown in Figure C.1 are the 
mass flowrates calculated by the reference standard at 0, 1, and 2GW/m’. There is 
surprisingly little difference between the simple model and the reference standard at these 


power densities. 
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Also note that the exit plenum O is adjusted from 2.66 to 1.10 as discussed in sec- 
tion 3.62. This gives a change in calculated flow (at 2 MW/m,) from 31 g/s (without the 
adjustment) to 37 g/s (=15% difference). The adjustment, however, did enable the 
behavior of the whole curve to be well represented. 

At this point, it must be noted that the savings in computer time is substantial. The 
time required for the reference standard to obtain the three data points shown on Figure 
C.1 is approximately 24 hours (8 hours per point) while the time required to generate the 
same data using the simple model is approximately 12 seconds. Further, the simple 


model is able to provide many more data points in the region of interest. 
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Figure C.1: Steady State Mass Flowrate vs Power Density 
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C.1.3 Power 
The rated output of the fuel element may be measured as the thermal power of the 
fuel element. This is calculated by multiplying the mass flowrate through the element by 


the change in enthalpy across the fuel element: 


Power = WAH Cal 


Figure C.2 shows the relationship between fuel element thermal power to power 
density and, as expected, the relationship is linear in nature and expresses a heat balance. 
The calculated value of thermal power at 2 MW/m’ is plotted for the reference model and 


corresponds to the value calculated by the simple model. 
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Figure C.2: Steady State Thermal Power vs Power Density 
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C.1.4 Exit Temperature 

Figure C.3 shows how the temperature of the coolant exiting the fuel element varies 
with power density. Because the exit plenum is modeled as a single control volume, the 
exit temperature represents an average temperature of the coolant exiting the element and 
does not take into account any temperature differences which may occur axially. The ref- 
erence standard has five axial subdivisions in the outlet plenum, each of which may have 
a different coolant temperature. Because of this axial separation and the ability of the 
reference standard to designate different axial fuel distributions, the temperatures calcu- 
lated in the exit plenum for the reference standard are higher at the lower axial position 
and decreases as the coolant exits the plenum. The simple model can not make this same 
representation and, as such, exit temperatures may be slightly different between the two 
models even though mass flowrates, power output and enthalpy changes may be the 
same. In this case, values for mass flowrate at various power density levels are near ref- 
erence standard values but are not exact. 

The reference temperature plotted in Figure C.3 is the average outlet plenum tem- 
perature at 2 and 2.1 MW/m’. This represents a difference of approximately 2.5% and 
should be considered acceptable for initial fuel element design but final analysis should 
be made by the more detailed model. 

Velocity is measured at the exit of the plenum. Because exit velocity of the coolant 
is inversely proportional to the density of the coolant, and hence directly proportional to 
the temperature of the coolant, it is not surprising to see the reference standard exit veloc- 
ity fall below the curve generated by the simple model. Exit velocity is shown as a func- 


tion of power density in Figure C.4. 
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Figure C.3: Steady State Exit Temperature vs Power Density 


90 





Steady State Exit Velocity 
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Figure C.4: Steady State Exit Velocity vs Power Density 
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C.1.5 Pressure 

The calculations of sections C.1 are based on constant inlet pressure and constant 
outlet pressure. As power increases, the total flow through the fuel element drops and the 
fraction of the total pressure drop increases in the control volumes at the higher tempera- 
tures (lower densities). C.5 shows these fraction by giving the calculated pressures at 
points located at the inlet and outlet of the element and at the inlet and outlet of the fuel. 
This also represents the boundaries of the three control volumes. 

As mentioned 1n T-1, the inlet plenum and the cold frit dominate the resistance at 0 
GW/m’ but as power density (and exit temperature) is increased, the outlet plenum 
assumes a greater share of the resistance. The resulting pressure drops across each con- 
trol volume are shown in Figure C.6. When totaled, the sum of the pressure drops will 
remain a constant 85 MPa, which is established as an initial condition. Figure C.6 also 
shows values for pressure drop across the outlet plenum in the reference standard. 
Although values at 0 and 1 GW/m’ are near to the values produced by the simple model, 
the value at 2 MW/m’ seems to be inconsistent. This inconsistency is unexplained and 1s 


judged to be unimportant when the consistent behavior of other variables 1s noted. 
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Figure C.5: Steady State Pressures vs Power Density 
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Figure C.6: Steady State Pressure Drops vs Power Density 
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C.2 NULL TRANSIENTS 

As a test for validating the simple model, null transients were analyzed at 0, 1, and 
2 MW/m’ and subsequent mass flowrate responses were examined. Figures C.7, C.8 and 
C.9 show that there 1s a smooth behavior throughout the transient and that there is negli- 
gible difference between the mass flowrate before the transient and the mass flowrate 
after the transient. Each of the transients start at .01 seconds and the scale for flowrate on 


the y-axis 1s expanded for better viewing. 
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Mass Flowrate Response to Null 
Transient at 0 GW/m3 for 
Baseline PIPE Fuel Element 


Flowrate (g/s) 
53.3 
53.2 
53.2 
53.1 
53.1 
0.01 0.012 0.014 0.016 0.018 0.02 


Time (s) 


Transient * Steady State 


Figure C.7: Mass Flowrate Response to Null Transient at 0 GW/m’* 





Mass Flowrate Response to Null 
Transient at 1 GW/m3 for 


Baseline PIPE Fuel Element 
Flowrate (g/s) 
44.40 
44.35 
44.30 
44.25 
44.20 
0.010 0.012 0.014 0.016 0.018 0.020 
Time (s) 


Transient x Steady State 


Figure C.8: Mass Flowrate Response to Null Transient at 1 GW/m’ 





Mass Flowrate Response to Null 
Transient at 2 GW/m3 for 


Baseline PIPE Fuel Element 
Flowrate (g/s) 
36.90 
36.85 
36.80 
36.75 
36.70 
0.010 0.012 0.014 0.016 0.018 0.020 
Time (s) 


Transient + Steady State 


Figure C.9: Mass Flowrate Response to Null Transient at 2 GW/m 





C.3 TRANSIENTS 

In Figures C.10, C.11, and C.12, calculations for a 1 second transient from 0 
GW/m’ to 2 GW/m’. Calculations for the simple model for mass flowrate, thermal 
power, and exit temperature are compared to the reference standard. As can be seen, the 
reference standard reaches the final set of values quicker than the simple model. Adyjust- 
ment could be made to correct the sluggishness of the simple model (e.g. by increasing 
the fuel-to-coolant heat transfer parameter and/or by artificially decreasing the mass of 
coolant in the fuel element). These adjustments have not been attempted in this thesis 
study. Both the simple model and the reference standard respond quickly; both end at 


essentially the same condition. 
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Mass Flowrate Response for Is 
Transient from 0 to GW/m3 for 
Baseline PIPE Fuel Element 
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Figure C.10: Mass Flowrate Response to Baseline 1 s Transient 
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Thermal Power Ouput Power for Is 
Transient from 0 to 2 GW/m3 for 
Baseline PIPE Fuel Element 


Power (kW) 


1000 


800 


600 


400 


200 





0.0 2.0 4.0 6.0 8.0 10.0 


Time (s) 


Nuclear Heat Gen Thermal Output (SM) 


* Thermal Output (RS) 





Transient Starts at t=1s 


Figure C.11: Thermal Power Response to Baseline 1 s Transient 
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Exit Temperature Response for 1s 
Transient from 0 to 2 GW/m3 for 
Baseline PIPE Fuel Element 


Temperature (K) 
2000 : 
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Transient Starts at t=ls 


Figure C.12: Exit Temperature Response to Baseline 1 s Transient 
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C.4 NUMERICAL STABILITY 

Accuracy of the solution of the control volume balance equations is illustrated in 
Figure C.13. That is, consider that the balance equations are written in terms of ordinary 
differential equations in time. Then consider solving these equations (as in the simple 
model) by discretizing the equations in the time variable. Figure C.13 gives the results 
related to the errors generated in this discretization. It does not address the question of 
discretization in space, however. 

The balance equations were discretized using timesteps of 100 ms, 10 ms, and 1 ms. 
All of the curves are very close and do not show much variation. The plot of mass flow- 
rate corresponding to timesteps of 1 ms and 10 ms seem to represent the solution of the 
differential equations better than the curve corresponding to 100 ms. The results for 100 


ms indicate an acceptable error (less than 5%) and could be used for faster calculations. 
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Mass Flowrate Response to Is 
Transient from 0 to 2 GW/m3 for 


Various Timesteps 





Flowrate (g/s) 
55.0 
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40.0 
35.0 
1.0 eS 2.0 2s 3.0 3.5 4.0 
Time (s) 


Timesteps 


dt=.01 s 





Transient Starts at t=1s 


Figure C.13: Mass Flowrate Response for Various Timesteps 
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APPENDIX D: 
SOURCE CODE AND DATA ENTRY 


D.1 SOURCE CODE DESIGN 

The source code for the analysis of the response of a particle bed reactor fuel ele- 
ment is provided in two parts: STEADY and UPPOWER. STEADY.C provides a means 
of calculating values at steady state conditions for a variety of fuel element 
configurations and operating parameters. This enables the user to analyze steady state 
conditions without having to run the transient program. UPPOWER.C is similar to 
STEADY.C but calculates transient values over a specified simulation period. 

Both source codes are written in the C programming language and compiled using 
Microsoft™ Quick-C. The programs are intended for use on a desktop personal computer 
equipped with a hard drive but may run on a computer equipped with floppy disk drives 
if one drive is designated as the "c drive". STEADY produces the output data file 
STEADYST.DAT and UPPOWER produces the data file POWER.DAT. The data files 
generated by STEADY and UPPOWER are unformatted files which are easily imported 
into a spreadsheet (ie, Lotus 123™) which then allows the user to review the output 
graphically. 

D.2 DATA ENTRY 

STEADY and UPPOWER provide fuel element geometry and initial conditions 

which are similar to the PIPE experiments as initial values of all program variables. 


These values may be changed via a series of input screens presented to the user prior to 
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beginning an analysis. The first screen (Figure E.1) initializes all important geometry 
variables which may be changed by responding to the prompts under the list of parame- 


ters. 
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THE FOLLOWING ARE PARAMETERS FOR 


Outer Radius 

Cold Frit Radius 

Cold Frit Thickness 

Dia of Cold Frit Particle 
Goid Frit: Porosity 
Particle Bed Thickness 


PPBR FUEL ELEMENT: 


.044580 m 
.029360 m 
.001870 m 
.00000270 m 
-685000 


012432 


mm 


-000500 
. 400000 
.000760 
. 230000 
. 950000 
. 888800 
. 265000 


Fuel Particle Diameter 
Particle Bed Void Fraction 
Hot Frit Thickness 

Hot Frit Porosity 

Inlet Manifold Factor 

Exit Manifold Factor 

Fuel Element Length 


™m 


own Ob WN b+ 


m 


oooococoococcocn$o 


= 
i 
= 
= 
= 


CHANGE A PARAMETER? (1=Y or O=N)1 


PLEASE ENTER THE PARAMETER NUMBER AND VALUE (eg 12, 
je aes er) 


Leto} 


CHANGE ANOTHER PARAMETER FOR ELEMENT A? ({1=Y or O=N)0 





Figure E.1: Data Input Screen #1 for STEADY 


Once the user is satisfied with the fuel element geometry, a listing of these values as 
they compare to the baseline PIPE fuel element 1s then shown on the screen (Figure E.2). 
This is followed by a screen which allows changes to specific fuel particle properties 
(Figure E.3) and a final input display which allows changes to general operating parame- 


ters (Figure E.4). 
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PARAMETER ELEMENT A BASELINE 


.044580 
.029360 
.001870 
.00000270 
-685000 
.012432 
-000500 
-400000 
.000760 
-230000 
- 950000 
- 888800 
-000000 m 


.044580 m 
.029360 m 
.001870 m 
-00000270 m 
.685000 
.012432 m 
.000500 m 
- 400000 
.000760 m 
.230000 

- 950000 

- 888800 
.265000 


Outer Radius 

Cold Frit Radius 

Cold Frit Thickness 

Dia of Cold Frit Particle 
Gold Frit Porosity 
Particle Bed thickness 
Fuel Particle Diameter 
Particle Bed Void Fraction 
Hot Frit Thickness 

fot Frit Porosity 

-Inlet Manifold Factor 
@2.Exit Manifold Factor 
13.Length of Element 


0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
0 
a 


Soro oqoonuqeec ee 


HIT ANY KEY TO CONTINUE 





Figure E.2: Data Input Screen #2 for STEADY 


Fuel Material 
Layer 1 Material 
Layer 2 Material 
Layer 3 Material 

Radius of Fuel 
Radius of Layer 1 
Radius of Layer 2 
Radius of Layer 3 

Density of Fuel 
10.Density of Layer i 
11.Density of Layer 2 
12.Density of Layer 3 
die Cp of Fuel 
14. of Layer 1 
5s of Layer 2 
16. of Layer 3 
17. k for Fuel 
18. k for Layer 1 
Lge k for Layer 2 
20. k for Layer 3 


Uc2 

Low D Car 

High D CarZrc 

yA Xe: 

0.000117 m 
0.000150 m 
0.000200 m 
0.000250 m 
10500.000000 kg/m3 
1000.000000 kg/m3 
1900.000000 kg/m3 
6300.000000 kg/m3 
200.000000 J/kgK 
3000.000000 J/kgK 
3000.000000 J/kgK 
200.000000 J/kgK 
30.000000 W/m2K 
1.500000 W/m2K 
3.000000 W/m2K 
40.000000 W/m2K 


So~yo ob Ww A 


oO 


og wnw vw Wd v # WH 8 8 Wh bo oP Pf oo 


CHANGE A PARAMETER? (1=Y or O=N)0 


CHANGE ANOTHER FARAMETER? (1=Y or O=N)0 





Figure E.3: Data Input Screen #3 for STEADY 
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THE FOLLOWING INITIAL AND FINAL CONDITIONS ARE SET 


a Initial Inlet Temperature = 300.000000 K 

2 Initial Inlet Pressure 2000.000000 KPa 
S. Initial Outlet Pressure 1915.000000 KPa 
4. Initial Power Density 0.000000 GW/m3 
5 
6 
ad 


Final Power Density 2.500000 GW/m3 
Power Density Increment 0.100000 sec 
Initial Guess at Mass Flowrate = 0.100000 kg/sec 
CHANGE A PARAMETER? (1=Y or O=N)1 


ENTER PARAMETER NUMBER, VALUE (eg 20, 40) 
3,1900 


CHANGE ANOTHER PARAMETER? (1=Y or O=N)0 





Figure E.4: Data Input Screen #4 for STEADY 


UPPOWER uses the same data input screens with the exception of input screen #4. 
UPPOWER allows the user to select a value for the initial and final value of inlet temper- 
ature, inlet pressure, outlet pressure and power density (Figure E.5). Further, the user 
may choose the duration of each transient. All transients start at the same time except for 
inlet temperature which has a 1 s delay to artificially simulate a preheating effect from a 


nuclear thermal rocket nozzle cooling jacket. 
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THE 


FOLLOWING INITIAL AND FINAL CONDITIONS ARE SET 


i. Initial Inlet Temperature = 300.000000 K 

ae Final Inlet Temperature = 300,000000 K 

3. Duration of Temperature Change = 1.000000 sec 

4. Initial Inlet Pressure = 2000.000000 KPa 
oe Final Inlet Pressure = 2000.,000000 KPa 
6. Duration of Pressure Change = 1.000000 sec 

wa Initial Outlet Pressure = 1915,000000 KPa 
B.. Final Outlet Pressure = 1915.000000 KPa 
re Duration of Pressure Change = 1.000000 sec 
EO. Initial Power Density = 0.000000 GW/m3 
ib es Final Power Density = 2.000000 GW/m3 
12.Duration of Power Density Change = 1.000000 sec 
Ss Time Delay to Transient = 1.000000 sec 
14. Time After Transient = 8,000000 sec 
LS. Time Step = 0,010000 sec 
aS. Initial Guess at Mass Flowrate = 0.100000 kg/sec 


CHANGE A PARAMETER? (1=Y or O=N)O 


CHANGE ANOTHER PARAMETER? (1=Y or O=N)0 
Data Sample Frequency? (print every xth data point...) 
10 


Figure E.5: Data Input Screen #4 for UPPOWER 


It should be noted that, as a final input, UPPOWER asks for a sampling frequency. 
This determines the amount of data written to the data file. Calculations are still per- 
formed at the required timestep, however. For example, if every data point were to be 
written to the data file using a timestep of 1 ms for a simulated transient lasting 10 s, 
there would be 10,000 data points written to the output file (one for each millisecond). If 
a sample frequency of 100 were to be used for the same transient, only 100 data points 
would be written to the output file (one every 100 ms). This serves only to reduce the 
data file to a manageable size and will not eliminate any of the calculations. 
D.3 Hydrogen Equations of State 

To analyze fuel element transients, properties of hydrogen coolant (such as density, 
viscosity and enthalpy) must be determined for various temperatures and pressures. 
Relationships to do this are placed at the end of the source code as subroutines and are 


accessed during the main part of the program as needed. 
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The density of hydrogen is determined by using a known value of hydrogen at 300 
K and 101.325 kPa (1 atmosphere) and using the ideal gas relationship to determine other 
values. The subroutine in the source codes determines the value of coolant density as a 
function of temperature and pressure. 

Viscosity, enthalpy, coolant heat conductivity and Prantdl number, on the other 
hand are determined as a function of temperature only since the contribution due to a 
change in pressure is negligible (R-1). To illustrate this, the enthalpy change from 200 K 
to 1200 K at a pressure of 101.325 kPa is 16406.2 kJ/kg and at a pressure of 2 MPa is 
16478.7 kJ/kg. The difference between the two values is approximately 0.4%. This is 
considered negligible for the sumple model. The other properties change in a similar 


fashion. 
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D.4 SOURCE CODE FOR STEADY 


/* STEADY.C */ 
/* Copywrite WILLIAM E. CASEY, MAY 1990. 


of 


/* The author hereby grants to MIT and to the US Government */ 
/* permission to reproduce and distribute copies of this code. */ 
/* This source code was written using the Microsoft QUICK-C  */ 


/* programming environment. =i 
/* Set up include files: */ 
#include <c:\msc\include\stdio.h> 

#include <c:‘\msc\include\conio.h> 

#include <c:\msc\include\graph.h> 

#include <c\msc\include\stdlib.h> 

#include <c‘\msc\include\math.h> 

#define PI 3.14159 

/* Define a structure for geometry related variables: 
struct fuelelem 


{ 


ss 


float OUTRAD; /* Outer radius of Inlet Plenum at 


float CFRAD; /* Quter radius of Cold Fnit 
float CFTHK; /* Cold Fnt Thickness 


* 


float CFDp; /* Diameter of Particles in Cold Frit */ 


float CFPOR; /* Cold Frit Porosity 
float PBEDTHK; /* Particle Bed Thickness 
floatPBEDDp; /* Diameter of Fuel Particles 


float PBEDVOID; /* Particle Bed Void Fraction 


float HFTHK; /* Hot Fnt Thickness 
float HFPOR; /* Hot Frit Porosity 


a 


a] 
sy 
ee 


oa 


os 


float MFIN; /* Manifold Factor forInlet Plenum */ 
float MFOUT; /* Manifold Factor for Outlet Plenum */ 


float LENGTH; /* Overall Length of Element * / 
i 
/* Define a structure for fuel particle composition: ia 
struct fuelpart 
char FUELMAT[10]; /* Material used for Fuel +) 
char LIMAT[10]; /* Material for Layer 1] =) 
char LIMAT[10]; /* Material for Layer 2 oy 
char L3MAT{[10]; /* Material for Layer 3 a7 
float FUELRAD; /* Radius of Fuel Core a 
float LIRAD; /* Radius of Layer 1 a) 
float L2ZRAD; /* Radius of Layer 2 */ 
float L3RAD; /* Radius of Layer 3 =) 
float FUELDEN; /* Density of Fuel | 
float L1 DEN; /* Density of Layer 1 oi 
float LZDEN; /* Density of Layer 2 | 
float L3DEN; /* Density of Layer 3 = 
float FUELCp; /* Cp for Fuel wl 
float L1 Cp; /* Cp for Layer 1 * / 
float L2Cp; /* Cp for Layer 2 7) 
float L3Cp; /* Cp for Layer 3 | 
float FUELK; /* Heat Transfer Coef for Fuel +) 
float LIK; /* Heat Transfer CoefforLayer1l */ 
float L2K; /* Heat Transfer CoefforLayer2 */ 
float L3K; /* Heat Transfer CoefforLayer3 */ 


lle 





/* Define a structure for initial and final conditions: ai! 
struct conditions 


floatINIT_TIN; /* Initial Value of T in */ 
float INIT_PIN; /* Initial Value of P in a) 
floatINIT_POUT; /* Initial Value of P out | 


float INIT_PRDEN; /* Initial Value of Power Density */ 
float FIN.PRDEN; /* Final Value of Power Density wi 
float delta_t; /* Increment a | 

float INIT_W; /* Initial Guess at Mass Flowrate */ 


e 

void showdefaults( struct fuelelem *e_ptr ); 

void fuelelemsum( struct fuelelem *e_ptr, struct fuelelem *f_ptr ); 
void showfuelpart ( struct fuelpart *g_ptr ); 

void showconditions ( struct conditions *h_ptr ); 


/* Define H2 state relationships found at end of source code: */ 
float h2density (float xx, float yy); 

float h2viscosity (float xx); 

float h2heatxfer (float xx); 

float h2cp (float xx); 

float h2prandle (float xx); 

float h2H (float xx); 


ies * BEGINNING OF COMPUTATIONAL CODE 2K 9k 2k 2K 9K 2k 9 2k 9 2K 9K 2K 2K 9k 2K ok 2K ok 2K 2k 2K 2k 2 ok aay 
main() 


/* Define variables for use in source code: sa) 
double WO, W1; 
double pinlI, pinl, pinIII, pout; 
double pbarl, pbarIl, pbarllI; 
double Tin, Tout, TbarlII; 
float PinI, PinII, PinIII, Pout, Pbarl, PbarlI1I, PbarII]I; 
float uin, uout, ubarlIl; 
float RI, RII, RIM, RTOTAL; 
float PDROPI1, PDROPII1, PDROPHI1; 
float PDROPI2, PDROPII2, PDROPIHII2; 
float HinI1, HinII1, HinlII1, Hout1; 
float HinI2, HinII2, HinlII2, Hout?2; 
float Hbarll, HbarI2, HbarII1, HbarlI2, HbarlII1, Hbarl2; 
float Qbed; 
float AinI, Aoutl, ACF, AinII, AoutI], AHF, Ain, AoutiIII; 
float AbarIl; 
float CPin, CPout, CPtemp; 
float Werror, CPerr; 
float Deq_in, Deq_out; 
float RIf, RIa, Rik, RIm; 
float RCF, Rel, Rell, RIlIa; 
float RUIf, RIMa, RIIk, RUIm; 
float mfout; 
float powerdensity, vout; 
int ch, chr; 
float aa, alpha, cc; 
float Htemp, Herror; 


Ls 





static struct fuelelem A = 

{ .04458, .02936, .00187, .0000027, .685, .012432, .000S, 
A, .00076, .23, .95, 1.1, .265 

Is 


Static struct fuelelem B = 
{ .04458, .02936, .00187, .0000027, .685, .012432, .0005, 
A, .00076, .23, .95, 1.1, .265 


I 

static struct fuelpart C = 

{ "“UC2", "Low D Car", "High D Car", "ZrC", .000117, .00015,\ 
.0002, .00025, 10500, 1000, 1900, 6300, 200, 3000, 3000, \ 
200, 30, 1.5, 3, 40 

Is 

static struct conditions D = 

{ 300, 2000, 1915, 0, 2.5, .1, .1 


Ie 

FILE *out; 

out = fopen ( "c:steadyst.dat","w+t" ); 
em= 1: 

while (ch == 1 ) 


showdefaults ( &A ); 
printf("“\nNCHANGE ANOTHER PARAMETER FOR ELEMENT A? (1=Y or 0=N)"); 
scanf("%i", &ch); 


} 

fuelelemsum ( &A, &B ); 

cn i: 

while (ch == 1 ) 

{ 
showfuelpart ( &C ); 
printf \“nNCHANGE ANOTHER PARAMETER? (1=Y or 0=N)"); 
scanf("%i", &ch); 

} 

eh = 1: 

while (ch == 1 ) 

{ 


showconditions ( &D ); 
printf("\nCHANGE ANOTHER PARAMETER? (1=Y or 0=N)"); 
scanf("%i", &ch); 

} 


/* Determination of Steady State Values-Initial Conditions */ 


fprintf( out,"STEADY STATE VALUES FOR THE FOLLOWING INITIAL CONDITIONS\o"); 
fprintf( out,"\nInlet Temp = %.2f K\n", D.INIT_TIN); 

fprintf( out,"Inlet Press = %.2f kPa\n", D.INIT_PIN); 

fprintf( out,"Outlet Press = %.2f kPa\n\n", D.INIT_POUT); 


fprintf(out, "PWRDN\Qbed\tW\NPDINPDINPDIINTouNNTbarINPinINPinINPinIIINtPoutNvout\n\n"); 
W1 =D.INIT_W; 
W0 = W1; 
Ain] = PI * (A.OUTRAD * A.OUTRAD - A.CFRAD * A.CFRAD); 
Aoutl = 2 * PI * A.CFRAD * A.LENGTH; 
ACF = AoutI * A.CFPOR; 
Ainl] = 2 * PI * (A.;CFRAD - A.CFTHK) * A.LENGTH * A.PBEDVOID; 
Aoutll = 2 * PI * (A.CFRAD - A.CFTHK - A.PBEDTHK) * A.LENGTH *\ 
A.PBEDVOID; 
AbarlIlI = .5 * (AinII + AoutII); 
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AinllI = 2 * PI * (A.;CFRAD - A.CFTHK - A.PBEDTHK - A.HFTHK) * \ 
A.LENGTH; 

AoutlII = PI * (A.CFRAD - A.CFTHK - A.PBEDTHK - A.HFTHK) *\ 
(A.CFRAD - A.CFTHK - A.PBEDTHK - A.HFTHK); 

AHF = AinlIIl * A.HFPOR; 

Deq_in = 4 * AinI /(2 * PI * (A-OUTRAD + A.CFRAD)); 

Deq_out = 4 * AoutIII / (2 * PI * (A-CFRAD - A.CFTHK - A.PBEDTHK -\ 
A.HFTHK)); 

for (powerdensity = D.INIT_PRDEN; powerdensity < D.FIN_PRDEN+.1;\ 

powerdensity += .1) 
{ 


Qbed = powerdensity * 1000000000 * PI * ((A.CFRAD - A.CFTHK) *\ 
(A.CFRAD - A.CFTHK) - (A.CFRAD - A.CFTHK - A.PBEDTHK) * \ 
(A.CFRAD - A.CFTHK - A.PBEDTHK)) * A.LENGTH; 

Tin = D.INIT_TIN; 

HinlI1 = h2H(Tin); 

Pinl = D.INIT_PIN; 

Pinll = D.INIT_PIN; 

Pinll] = D.INIT_POUT; 

Pout = D.INIT_POUT; 

Tout = Tin; 


WO0 =WI1; 
Houtl = Qbed/W0O + Hinl1; 
i 
Htemp = h2H(Tout); 
Herror = Houtl - Htemp; 
Tout = Tout + Herror/100000; 
} while ( Herror/1000000 > .000001 ); 
i Set up all variables yf! 
TbarlII = .5*(Tin + Tout); 
Pbarl = .5*(PinI + Pinll); 
Pbarll = .5*(PinII + PinII); 
Pbarlll = .5*(PinIlI + Pout); 
pbarl = h2density(Tin, Pbar!); 
pbarlII = h2density(TbarIlI, PbarlII); 
pbarIII= h2density(Tout, PbarIIT); 
uin = h2viscosity(Tin); 
ubarll = h2viscosity(Tbarl]); 
uout = h2viscosity(Tout); 
pinll = h2density(Tin, PinI); 
pinllI = h2density(Tout, PinIII); 
pout = h2density(Tout,Pout); 
[is Inlet Resistance Calculations ) 
Rel = W0*Deq_in / (uin* Ainl); 
RIf = .138*pow(Rel, -.151)*(A.LENGTH/Deq_in)/((2*pbarl* AinI* Ainl); 
Rla = (1/(AinII*AinII) - 1/(AinI* AinI))/(2*pbarl); 
RIk = (1-(AinI/AoutI))*(1-(AinI/A outl))((2*pbarI*AinI*AinI) \ 
+ .5 * (1-(ACF/Aoutl))*(1-(ACF/Aoutl)) / (2*pinII*ACF*ACF) \ 
+ .5 * (1-(AinII/ACF))*(1-(AinII/ACF)) / (2*pinII* AinII* AinID); 
RIm = A.MFIN / (pbarI* AinI* AinI); 
RIm = RIm + A.MFIN/(pbarI* AoutI*Aout!); 
RCF = 150*uin*(1-A.CFPOR)*(1-A.CFPOR)/(A.CFDp*A.CFDp*pinlI* AoutI* W0* \ 
A.CFPOR*A.CFPOR*A.CFPOR) \ 


|B 





- TT 





“s : j 


+ 1.75*(1-A.CFPOR)/(A.CFDp*pinlIl*Aout!* Aoutl*A.CFPOR*A.CFPOR* \ 
A.CFPOR); 
Ree = RCFE * A.CFTHK: 
RI = RIf/2 + RIla + RIk + Rim + RCF; 
/* Resistance Calculations for the Particle Bed 3 
RII = 150*ubarlII*(1-A.PBEDVOID)*(1-A.PBEDVOID)/(A.PBEDDp*A.PBEDDp* \ 
pbarlII* AbarlII* WO*A.PBEDVOID*A.PBEDVOID) wt 
+ 1.75*(1-A.PBEDVOID)/(A.PBEDDp*pbarlII*AbarlI*AbarlI*A.PBEDVOID); 
RII = RIT * A,PPBEDTHK; 
Rila = ( 1/(pinlII* AinlI) - 1/(pinI]*AinI) )*((AinII+AinIID)/ \ 
(2* AinII* AinllI)); 
RII = RII + Ria; 
[* Resistance for Exit Plenum | 
RellI = W0*Deq_out / (uout* AoutlII); 
RIIIf = .138*pow(Relll, -.151)*(A.LENGTH/Deq_out)/(2*pbarlII*AoutlII \ 
* Aout); 
Rilla = (1/(AoutII*A out T)-1/(AoutI* Aoutl!))/(2*pbarlIlI); 
RITIk = .5*(1-AHF/AoutII)*(1-AHF/A outIT)/(2*pinllI* AHF*AHF) \ 
+ (1-AHF/AinIIT)*(1-AHF/AinIIT)/(2*pinlI*AHF*AHF) \ 
+ (1-AoutITI/AinII)*(1-AoutITI/Aini)/(2*pbarlII* \ 
AoutIII* AoutIII); 
mfout = AMMFOUT; 
RiIm = mfout / (pbarIII*A outI*AoutIII) + mfout / (pbarlII\ 


* Ainil* AinlII); 
RITI = RUOIf/2 + RUla + RIUIk + RITIm; 
= Final Assembly aa 


RTOTAL = RI + RH + RII; 
W1 =(PinI - Pout)*1000/RTOTAL; 
W1 = sqrt(W1); 
Werror = WO - W1; 
Werror = fabs( Werror); 
PDROPI1 = RI*W1*W1/1000; 
PDROPII1 = RII*W1*W1/1000; 
PDROPHI1 = ROI*W1*W1/1000; 
Pin = PinI-PDROPI1; 
PinI] = Pout+PDROPIT1; 
vout = W1/(pout*AoutIID); 

} while (Werror > .0000001); 
W1=W1*1000; 
Qbed = Qbed/1000; 


fprintf(out,"%.2N%.1N%.SN%.2N%.2N%G.2NGD.2NG.2NG.2NG.2NG.2NG.2N%G.2fn", \ 
powerdensity, Qbed, W1,PDROPI1,PDROPI1,PDROPIII1 , Tout, TbarlII, Pi- 
nI,PinII,PinlIII,Pout,vout); 
printf ("Pwrden = %Z{NW = %f\n"", powerdensity, W1); 
Qbed = Qbed* 1000; 
W1=W1/1000; 
if (Tout > 2700) 
goto meltdown; 
} 
meltdown: 
if (Tout>2700) 
{ printf("Tout > 2700"); 
fprintf(out,"Temp exceeds 2700"); 
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Hie Be OK 2h OK 2k 2k 2k 2K 2k 2 2k 2k 2k 2k 2k 2k 2k 2k 2k 2k 2k 26 2k 3K 2K 26 25 2 2k ok 2 26 26 3K 3 26 2K 2 ok 2k 2K 2 26 2k 2k 2 26 2k 2k 2 2 26 2k 2 2k 2 2 26 3 26 2k 2 2 2c ok 2k “y 
void showdefaults( struct fuelelem *e_ptr ) 


{ 


int num; 

float entry; 

mt c: 

_clearscreen(Q); 

prntf("THE FOLLOWING ARE PARAMETERS FOR PPBR FUEL ELEMENT "n\n"); 
printf(''1. Outer Radius = %f m\n", e_ptr->OUTRAD); 

printf("2. Cold Frit Radius = %f m\n", e_ptr->CFRAD); 

printf(''3. Cold Frit Thickness = %f m\n", e_ptr->CFTHK); 

prntf("4. Dia of Cold Frit Particle = %.8f m\n", e_ptr->CFDp); 

printf("'S. Cold Frit Porosity = %f\n"", e_ptr->CFPOR); 


printf(''6. Particle Bed Thickness = %f m\n", e_ptr->PBEDTHK); 
printf(''7. Fuel Particle Diameter = %f m\n", e_ptr->PBEDDp); 
printf("8. Particle Bed Void Fraction = %f\n", e_ptr->PBEDVOID); 
printf('"9. Hot Frit Thickness = %f m\n", e_ptr->HFTHK); 
printf("10. Hot Frit Porosity = %f\n"", e_ptr->HFPOR); 
printf("11. Inlet Manifold Factor = %f\n", e_ptr->MFIN); 
printf("12. Exit Manifold Factor = %f\n", e_ptr->MFOUT); 
printf("13. Fuel Element Length = %f m\n", e_ptr->LENGTH); 
printf(‘\n\n\nCHANGE A PARAMETER? (1=Y or 0=N)"); 
scanf( "%i", &c ); 
if(c==0) 

returm; 
printf(\nPLEASE ENTER THE PARAMETER NUMBER AND VALUE (eg 12, 1.705)Nn"); 
scanf( "%i, %f", &num, &entry ); 
switch (num) 


case 1: 
e_ptr->OUTRAD = entry; 
break; 

ease 2: 
e_ptr->CFRAD = entry; 
break; 

case 3: 
e_ptr->CFTHK = entry; 
break; 

case 4: 
e_ptr->CFDp = entry; 
break; 

case 5: 
e_ptr->CFPOR = entry; 
break; 

case 6: 
e_ptr->PBEDTHK = entry; 
break; 

case 7: 
e_ptr->PBEDDp = entry; 
break; 

case 8: 
e_ptr->PBEDV OID = entry; 
break; 

case 9: 
e_ptr->HFTHK = entry; 
break; 

case 10: 


Diy 





e_ptr->HFPOR = entry; 
break; 

case ll: 
e_ptr->MFIN = entry; 
break; 

case 12: 
e_ptr->MFOUT = entry; 
break; 

Gase 13: 
e_ptr->LENGTH = entry; 
break; 

default: 
printf(\n\nTRY ANOTHER PARAMETER\n"); 
_clearscreen(0); 
break; 

} 


retum,; 


ig OR 2K 2K 2k 6 2k 2 3k 2k 2k 2K 2k 2 3 2K ok 26 2K 2k 2k 2 2K 2k 2 ok 2K 2k 2 oo 2K 2 2K 2 ok 2 2k 26 2K ok ok oi 2K 2K ok OK ok ok 2k ok 2K 2 2 ok 2K 2 2K ok ok 2k 3K 2k 2K ok ok ok a 
void fuelelemsum ( struct fuelelem *e_ptr, struct fuelelem *f_ptr ) 


{ 
_clearscreen (0); 
printf ("PARAMETER\NNELEMENT A\NBASELINE\n\n"); 
printf("1. Outer Radius WN%AN%E m\n", e_ptr->OUTRAD, f_ptr->OUTRAD); 
printf("2. Cold Frit Radius\W%fANNt%f m\n", e_ptr->CFRAD, f_ptr->CFRAD); 
printf("3. Cold Frit Thickness\N%ANI%f m\n", e_ptr->CFTHK, f_ptr->CFTHK); 
printf("4. Dia of Cold Frit ParticleN%.8NNN%.8f m\n", e_ptr->CFDp, f_ptr->CFDp); 
pmntf("5. Cold Frit PorosityWN%fN\Nt%fin"", e_ptr->CFPOR, f_ptr->CFPOR); 
printf("6. Particle Bed thickness %f\N%f m\n", e_ptr->PBEDTHK, f_ptr->PBEDTHK); 
printf("7. Fuel Particle DiametertZfN\N%f m\n", e_ptr->PBEDDp, f_ptr->PBEDDp); 
printf('"8. Particle Bed Void FractionNZf\N%fn", e_ptr->PBEDVOID, f_ptr->PBEDVOID); 
printf("9. Hot Frit Thickness\WN%f\N%f m\n", e_ptr->HFTHK, f_ptr->HFTHK); 
printf("10.Hot Frit Porosity\N%f\N% fn", e_ptr->HFPOR, f_ptr->HFPOR); 
pnintf("11.Inlet Manifold FactoNw%fhN%fn", e_ptr->MFIN, f_ptr->MFIN); 
printf("12.Exit Manifold FactoWWN@AN%fn",", e_ptr->MFOUT, f_ptr->MFOUT); 
printf("13.Length of ElementNN%f m\N%f\n"," e_ptr->LENGTH, f_ptr->LENGTH); 
printf(‘\nnHIT ANY KEY TO CONTINUE"); 
getch(); 


Io 3K 2K 2K 2 2K 2K 2 2 2k 2 2k 2k 2 2k 26 2K og 2K 2k 3 2k 2k 2k 2k 2k ok ok 2k 2k 2 2k 2 2K 2 2K 2k 2k ok ok 26 oi 26 2g 2k 2k 2k 2k 2k 2k 2k 2K ok 2 2k ok 2k 2K 2k 2k ok 2k 2k 2 2k 2K 2k * / 
void showfuelpart ( struct fuelpart *g_ptr ) 


int num; 

float entry; 

int C; 

_clearscreen(Q); 

printf("THE FOLLOWING ARE PARAMETERS FOR THE FUEL PARTICLE\n"); 
printf("l. Fuel Material = %s\n", g_ptr->FUELMAT); 
printf("2. Layer 1 Material = %s\n", g_ptr->L1MAT); 
printf("3. Layer 2 Material = %s\n", g_ptr->L2MAT); 
printf("4. Layer 3 Material = %s\n", g_ptr->L3MAT)); 
pnntf("5. Radius of Fuel = %f m\n", g_ptr->FUELRAD); 
printf("6. Radius of Layer 1 = %f m\n", g_ptr->L1 RAD); 
printf("7. Radius of Layer 2 = %f m\n", g_ptr->L2RAD); 
printf("8. Radius of Layer 3 = %f m\n", g_ptr->L3RAD); 
printf("9. Density of Fuel = %f kg/m3\n", g_ptr->FUELDEN); 
printf("10.Density of Layer 1 = %f kg/m3\n", g_ptr->L1 DEN); 
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printf("11.Density of Layer 2 = %fkg/m3\n", g_ptr->L2DEN); 
printf("12.Density of Layer 3 = %f kg/m3\n", g_ptr->L3DEN); 
printf("13. Cp of Fuel = %f J/kgK\n"," g_ptr->FUELCp); 
printf("14. Cp of Layer 1 = %f J/kgK\n", g_ptr->L1Cp); 
printf("15. Cp of Layer 2 = %f J/kgK\n", g_ptr->L2Cp); 
printf("16. Cp of Layer 3 = %f J/kgK\n", g_ptr->L3Cp); 
printf("17. k for Fuel = %f W/m2K\n", g_ptr->FUELK); 
printf('18. k for Layer 1 = %f W/m2K\n", g_ptr->L1K); 
printf("19.  k for Layer 2 = %f W/m2K\n", g_ptr->L2K); 
printf("20. k for Layer 3 = %f W/m2K\n", g_ptr->L3K); 
printf(\nCHANGE A PARAMETER? (1 =Y or 0=N)"); 

scanf( "%i", &c ); 

nm (c=—0 ) 


retum,; 
printf(\nENTER PARAMETER NUMBER ( NOT 1-4 ), VALUE (eg 20, 40)\n"); 
scanf( "%1, %f", &num, &entry ); 
switch (num) 


case 5: 
g_ptr->FUELRAD = entry; 
break; 

case 6: 
g_ptr->L1RAD = entry; 
break; 

Case 7. 
g_ptr->L2RAD = entry; 
break; 

case 8: 
g_ptr->L3RAD = entry; 
break; 

case 9: 
g_ptr->FUELDEN = entry; 
break; 

case 10: 
g_ptr->L1DEN = entry; 
break; 

case 11: 
g_ptr->L2DEN = entry; 
break; 

case 12: 
g_ptr->L3DEN = entry; 
break; 

case 13: 
g_ptr->FUELCp = entry; 
break; 

case 14: 
g_ptr->L1Cp = entry; 
break; 

case 15: 
g_ptr->L2Cp = entry; 
break; 

case 16: 
g_ptr->L3Cp = entry; 
break; 

casera. 
g_ptr->FUELK = entry; 
break; 

case 18: 
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g_ ptr->LIK = entry; 
break; 

case 19: 
g ptr->L2K = entry; 
break; 

case 20: 
g_ ptr->L3K = entry; 
break; 

default: 
prntf(“CANNOT CHANGE THIS PARAMETER...SORRY.\n"); 
getch(); 
break; 

} 


retum; 


[hg 3K 2k 2k 2k 2k 2k 2k 2k 2k 2k 2k 2k ok ok ok ok ok 2c ok ok ok 2k ok 2k kc ok ak 2k ok ok ok 2k 2k 2c ok 2k ok 2k akc ok akc ok ok ok ak ok ok 2k ok 2k ac ok 2k akc ok 2k 2k ok 2k akc ok ok akc ok ak 2k -/ 
void showconditions ( struct conditions *h_ptr ) 


int num; 

float entry; 

int C; 

_clearscreen(0); 

printf("THE FOLLOWING INITIAL AND FINAL CONDITIONS ARE SET \n\n"); 
printf("1. Initial Inlet Temperature = %f K\n", h_ptr->INIT_TIN); 


printf("'2. Initial Inlet Pressure = %f KPa\n", h_ptr->INIT_PIN); 
printf("3. Initial Outlet Pressure = %f KPa\n", h_ptr->INIT_POUT); 
printf("'4. Initial Power Density = %f GW/m3\n", h_ptr->INIT_PRDEN); 
printf('S. Final Power Density = %f GW/m3Nn", h_ptr->FIN_PRDEN); 


printf("6. Power Density Increment = %f sec\n", h_ptr->delta_t); 
printf(""7. Initial Guess at Mass Flowrate = %f kg/sec\n", h_ptr->INIT_W); 
printf \nCHANGE A PARAMETER? (1=Y or 0=N)"); 
scanf( "%i", &c ); 
if (c==0) 

retum; 
printf(\nENTER PARAMETER NUMBER, VALUE (eg 20, 40)\n"); 
scanf( "%i, %f", &num, &entry ); 
switch (num) 


case 1: 
h_ptr->INIT_TIN = entry; 
break; 

case 2: 
h_ptr->INIT_PIN = entry; 
break; 

case 3: 
h_ptr->INIT_POUT = entry; 
break; 

case 4: 
h_ptr->INIT_PRDEN = entry; 
break; 

case 5: 
h_ptr->FIN_PRDEN = entry; 
break; 

case 6: 
h_ptr->delta_t = entry; 
break; 

case 7: 
h_ptr->INIT_W = entry; 
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break; 
default: 
pnntf("CANNOT CHANGE THIS PARAMETER...SORRY An"); 
getch(); 
break; 
} 


retum, 


1 2K 2K 2 2 2k oie 2k oie 2k ak ic 2k 2K 6 go ic 6 2k IK 2k 2k 2 2k 2k 26 IK 2k 2k 2c 2k ak 2k oi ic ak 2k 2k 2k 2k ik 2k 2c 2k 2k 2k 2k 2k oie 2c 2 2k 2k 2k aie 2k 2s ic ak ic ak 2s 2s 2 of oc 7 


float h2density(float xx, float yy) 


float dens; 
dens = .08185 * 300 / xx * yy / 101.325; 
retum dens; 


/* 2K 2 2 2k 2k 2k 2k 2k ik 2K 2k 2k 2K 2 2k ok 2k 2k 2K 2 2k 2k 2k 2k 2k ak 2K 2k 2k 2k 2k 2k 2k 2k 2k 2c 2k 2k 2 oie 2K 2c 2k 2K 2k 2 2k 2K 2k 2K 2K 2 2k 2 2k 2k 2c 2K 2 2k 2k 2k 2k 2k 2k 2k | 
float h2viscosity(float xx) 


float visc; 


visc = .000008963 * pow( xx/300, .6733); 
return visc; 


/* 2K OK 2 oie 2 2h 2k 2k 2k 2k 2k oie 2c 2 2c 2 2c 2k oie ik ak 2k 2k oie ak ak 2c 2k 2k 2k oc ofc ok okt 2k 2c 2k 2c 2k oie ik 2c 2k ok ofc oe oc ai aK ok aK 2k 2k oi oie ofc oc 2k oie aK oe ofc kc 2k ak ok aE 
float h2heatxfer(float xx) 


float tem1; 

int 1; 

float tem2; 

float heat[45] = {0,.0362, .0665, .0981, .1282, .1561, .182, .206, .228\ 
le 272.292, 315, 5339-551, 5669, 384, .398,\ 
412, .426, .44, .452, .464, .476, .488, .500, 512, \ 
524, .536, 548, 560, .572, .584, .596, .608, .62, \ 
.632, .644, .656, .668, .680, .692, .704, .716, .728 }; 

tem2 = xx / 50.0; 

i = floor(tem2); 
tem1 = heat[i] + (heat[i+1] - heat[i]) * ((xx - Gi * 50)) /50); 
return tem1; 


/* 3 AK 2 AK 2k 2k 2k 2k 2c 2c 2 oie 2c 2k 2 ok 2c 2 2c ik akc ic oc kc ic 2 2 2k ofc 2k aie 2c 2k 2k 2c 2k 2k ic kc ic afk oie aie ic 2c 2k aie aie 2k 2c oie oie ofc kc is 2k oie 2k 2k ic oc 2k aK 2 ok * / 


float h2cp(float xx) 
{ 


float temcp; 
if (xx <= 420 ) 
temcp = 4.1868 * 1000 * (1.5395 + .0150825 * xx - 4.02449e-5 *\ 
xx * xx + 3.63544e-8 * xx * xx * xx); 
else 
temcp = 4.1868 * 1000 * (3.58927 - 5.55096e-4 * xx + 6.94235e-7 * xx\ 
* xx - 1.45155e-10 * xx * xx * xx); 
return temcp; 


/* 9K 2k 2k 2K 2K 2 2K 2k 2K 2 2K A 2k 2K 2k 2K 2c 2 2 ok 2K 2 2k 2k 2k 2 2 2 2k 2k 2k 2 2 ok 2k ok 2k 2k 2K ok 2k 2 2K 2k 2k 2 2k 2K 2k 2k 2 2 2 2K 2c 2K ok 2K 2k ok 2k 2k 2k 2K 2k ok ll 


float h2H (float xx) 
{ 


float temH; 
float temH1; 
temH=0; 
temH 1=0; 


| 





if (xx <= 420) 


{ 

temH=4.1868*1000*(1.5395*xx + .0150825*xx*xx/2 - 4.02449e-5* x 
XX*xx*xx/3 + 3.63544e-8*xx *xXx*xXxX*xx/4); 

temH=temH - 4.1868*1000*(1.5395*50 + .0150825*50*50/2 - 4.02449e-5* \ 
50*50*50/3 + 3.63544e-8*50*50*50*50/4); 

} 


else 


{ 

temH1=4.1868*1000*(3.58927*xx - 5.55096e-4*xx*xx/2 + 6.94235e-7*  \ 
XX *xx*xx/3 - 1.45155e-10*xx*xx*xx*xx/4); 

temHl=temH1 - 4.1868*1000*(3.58927*420 - 5.55096e-4*420*420/2+ \ 
6.9423 5e-7*420*420*420/3 - 1.4515Se-10*420*420*420*420/4); 

temH=temH1 + 4.1868*1000*(1.5395*420+.0150825*420*420/2-4.02449e-5* \ 
420*420*420/3 + 3.63544e-8*420*420*420*420/4); 

temH=temH - 4.1868*1000*(1.5395*50 + .0150825*50*50/2 - 4.02449e-5* \ 
50*50*50/3 + 3.63544e-8*50*50*50*50/4); 


return temH; 


iba 2K 2 2K 2 2k 2k 2k 2 2k 2k 2k 2 2k 2K 2k 2k 9 2K 2 2K 2 2k 2k 2K oie 2g 2k 2k 2 2 3k 2K 2K 2k 2k 2k 2 2k 2K 2k 2k 2k 2k 26 2K 2 2 2K 2 2K 2 2K 2 2K 2 26 2 2K 2K 2K 2K 2K 2K 2K 3 2k */ 


float h2prandle(float xx) 
{ 


float tem1; 

int 1; 

float tem2; 

float pran[23]= {.000, .712, .719, .706, .690, .675, .664, \ 
.659, .664, .676, .686, .703, .715, .733, \ 
.748, .763, .778, .793, .808, .823, .838, \ 
.853, .868}; 

tem2 = xx / 100.0; 

1 = floor(tem2); 
tem1 = pran{i] + (pran[i+1] - pran[i]) * ((xx - (i * 100)) / 100 ); 
return tem1; 


‘haa 2K 2K 2 2k 2 2 2 2k 2 ok 2 2 2 2 2 2k 2k 2 2k 2k ok 2k 2k 2 2k 2 2k 2K 2 2k 2k 2 2K 2K 2K 2K 2 2k 2k 2k 2k 2 2 a 2k 2 2 2 2k 2k 2 kk 2K 2 ok ok 2 2k 2K ok 2 * / 
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E.5 SOURCE CODE FOR UPPOWER 


/* UPPOWER.C 2 

/* Copywrite WILLIAM E. CASEY, MAY 1990. a 
/* The author hereby grants to MIT and to the US Government */ 
/* permission to reproduce and distribute copies of this code. */ 

/* This source code was written using the Microsoft QUICK-C  */ 
/* programming environment. ay 

/* Set up include files: a 

#include <c:‘\msc\include\stdio.h> 

#include <c:\msc\include\conio.h> 

#include <c:\msc\include\graph.h> 

#include <c:\msc\include\stdlib.h> 

#include <c\msc\include\math.h> 

#define PI 3.14159265442 


struct fuelelem 
float OUTRAD; /* Outer radius of Inlet Plenum * / 
float CFRAD; /* Outer radius of Cold Frit ai 
float CFTHK; /* Cold Frit Thickness *) 
float CFDp; /* Diameter of Particles in Cold Frit */ 
float CFPOR; /* Cold Frit Porosity a 
floatPBEDTHK; /* Particle Bed Thickness “| 


float PBEDDp; /* Diameter of Fuel Particles | 
float PBEDVOID; /* Particle Bed Void Fraction ag 
float HFTHK; /* Hot Frit Thickness * | 

float HFPOR; /* Hot Frit Porosity */ 

float MFIN; /* Manifold Factor for Inlet Plenum */ 
float MFOUT; /* Manifold Factor for Outlet Plenum */ 


float LENGTH; /* Overall Length of Element +} 
}; 
struct fuelpart 
char FUELMAT[10]; /* Material used for Fuel oy 
char LIMAT[10}; /* Material for Layer 1 st) 
char LZMAT[10]; /* Material for Layer 2 =) | 
char L3MAT[10}; /* Material for Layer 3 i) 
float FUELRAD; /* Radius of Fuel Core a, 
float LIRAD; /* Radius of Layer 1 a 
float LZARAD; /* Radius of Layer 2 =o)! 
float L3RAD; /* Radius of Layer 3 all 
float FUELDEN; /* Density of Fuel t/ 
float LIDEN; /* Density of Layer 1 / 
float L2DEN; /* Density of Layer 2 i 
float L3DEN; /* Density of Layer 3 oy) 
float FUELCp; /* Cp for Fuel 7 
float L1Cp; /* Cp for Layer 1 ny 
float L2Cp; /* Cp for Layer 2 i) 
float L3Cp; /* Cp for Layer 3 sa | 
float FUELK; /* Heat Transfer Coef for Fuel Sf 
float L1K; /* Heat Transfer Coef for Layer1 */ 
float L2K; /* Heat Transfer Coef forLayer2 */ 
float L3K; /* Heat Transfer CoefforLayer3 */ 


Ie 


struct conditions 
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float INIT_TIN; /* Initial Value of T in a 


float FIN_TIN; /* Final Value of T in | 
float DUR_TIN;  /* Duration of T in Transient “/ 
float INIT_PIN; /* Initial Value of P in * / 
float FIN PIN;  /* Final Value of P in | 
float DUR_PIN; /* Duration of P in Transient “ay | 
float INIT_POUT; /* Initial Value of P out | 
float FIN-POUT; /* Final Value of P out x) 


float DUR_POUT; /* Duration of P out Transient */ 

float INIT_PRDEN; /* Initial Value of PowerDensity */ 
float FIN_PRDEN; /* Final Value of Power Density */ 
float DUR_PRDEN; /* Duration of Power Density Transient */ 
floatt_DELAY; /* Time Delay prior to Transient wl 
floatt_AFTER; /* Window of Time After Transient */ 
float delta_t; /* Time Step a 

float INIT_W; /* Initial Guess at Mass Flowrate */ 


3 

void showdefaults( struct fuelelem *e_ptr ); 

void fuelelemsum( struct fuelelem *e_ptr, struct fuelelem *f_ptr ); 
void showfuelpart ( struct fuelpart *g_ptr ); 

void showconditions ( struct conditions *h_ptr ); 

double h2density(double xx, double yy); 

double h2viscosity(double xx); 

double h2heatxfer(double xx); 

double h2cp(double xx); 

float h2prandie(float xx); 

float h2H(float xx); 

i* 2K 3c BEGINNING OF COMPUTATIONAL CODE 2K ok 2 2 2K 26 ok ok ok ak a ak ok ok ok ac 2k ak ak “a 
main() 


/* Define variables: af 
int ch, chr; 
float W1, W2; 
float W; 
double Werror; 
float pinl, pinlI, pinlII, pout; 
float pbarlI, pbarlII, pbarIII; 
float uin, uout, ubarII; 
float Tin, TinI, TinII, TinIII, Tout, TbarlI, TbarII, TbarlIII; 
float PinI, PinII, PinIII, Pout, Pbarl, PbarlII, PbarIII; 
float Pinl1, PinII1, PinIII1, Poutl, PbarI1, Pbarl{1, PbarIII1; 
float PDROPI1, PDROPII1, PDROPIII1, PDROPI2, PDROPII2, PDROPIII2; 
float Hinl, HinI1, HinlI1, HinlII1, Hout1; 
float Hin2, HinI2, HinII2, HinIII2, Hout2, Htemp; 
double Herror; 
float Deq_in, Deq_out; 
float AinI, AoutI, ACF, AinII, AoutII, AHF, AinIII, AoutIII; 
float AbarlI; 
float voll, voll], vollI; 
float RI, RII, RIM, RTOTAL; 
float RIf, RIla, RIk, RIm, RIa; 
float RCF, Rel, Rell, Rell]; 
float RIJIf, RIIIa, RIJk, RIIIm, mfout; 
float Qbed, powerdensity, vout, power, powerout; 
double time, transtime, deltat; 
float mal, ma2, ma3, ma4, mat; 
float U1, U2, U3, U4, UT, UBAR, Cpbar; 
float f1, f2, £8, fbar; 


124 





float Nu, Pr, k, h; 

float I1, 12, 13, M1, M2, M3; 

float AA, BB, DD, EE, N; 

float Vcv, Av; 

float Tf2, Tf1, Qs; 

float WA,pfuelave; 

static struct fuelelem A = 

{ .04458, .02936, .00187, .0000027, .685, .012432, .000S, 
4, .00076, .23, .95, 1.1, .265 

b 


static struct fuelelem B = 

{ .04458, .02936, .00187, .0000027, .685, .012432, .000S, 
4, .00076, .23, .95, 1.1, .265 

i 


static struct fuelpart C = 

{ "UC2", "Low D Car", "High D Car", "ZrC", .000117, .00015,\ 
.0002, .00025, 10500, 1000, 1900, 6300, 200, 3000, 3000, \ 
200, 30, 1.5, 3, 40 

i 

static struct conditions D = 

{ 300, 300, 1, 2000, 2000, 1, 1915, 1915, 1, 0, 2, 1, 1, 8, .01, .1 

FILE *out; 

out = fopen ( "c:power.dat","w+"' ); 


Cni= |: 
while (ch = 1 ) 
{ 


showdefaults ( &A ); 
printf \nCHANGE ANOTHER PARAMETER FOR ELEMENT A? (1=Y or 0=N)"); 
scanf("%1", &ch); 


} 

fuelelemsum ( &A, &B ); 
cu — |: 

while (ch == 1 ) 


showfuelpart ( &C ); 
printf’ nCHANGE ANOTHER PARAMETER? (1=Y or 0=N)"); 
scanf("%i", &ch); 
} 

en = ic 

while (ch == 1 ) 

{ 
showconditions ( &D ); 
printf’ nCHANGE ANOTHER PARAMETER? (1=Y or 0=N)"); 
scanf("%i", &ch); 


} 

printf("Data Sample Frequency?(print every xth data point..)\n"); 

scanf("%f", & EE); 

fprintf( out,"VALUES FOR THE FOLLOWING INITIAL CONDITIONS‘n"); 
fprintf( out,"\n_ Initial Inlet Temp = %.2f K\n", D.INIT_TIN); 

fprintf( out,” Initial Inlet Press = %.2f kPa\n", D.INIT_PIN); 

fprintf( out,” Initial Outlet Press = %.2f kPa\n\n", D.INIT_POUT); 

fprintf(out, "time\ttr MWHI2NNHIT2\NHITI2\NtHout2\tpden\tW2\n"); 
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W1=D.INIT_W;: 

Ain] = PI * (A.OUTRAD * A.OUTRAD - A.CFRAD * A.CFRAD); 

Aout] = 2 * PI * ACCFRAD * A.LENGTH; 

ACF = AoutI * A.CFPOR;: 

Ainl] = 2 * PI * (A.CFRAD - A.CFTHK) * A.LENGTH * A.PBEDVOID; 

AoutlI = 2 * PI * (A.CFRAD - A.CFTHK - A.PBEDTHK) * A.LENGTH *\ 

A.PBEDVOID; 
AbarlIl = .5 * (AinII + AoutlI); 
AinllI =2 * PI * (A.CFRAD - A.CFTHK - A.PBEDTHK - A.HFTHK) *\ 
A.LENGTH; 

AoutlII = PI * (A.CFRAD - A.CFTHK - A.PBEDTHK - A.HFTHK) *\ 
(A.CFRAD - A.CFTHK - A.PBEDTHK - A.HFTHK); 

AHF = AinllI * A.HFPOR; 

Deq_in = 4 * Ain] /(2 * PI * (A.OUTRAD + A.CFRAD)); 

Deq_out = 4 * AoutIII / (2 * PI * (A.CFRAD - A.CFTHK - A.PBEDTHK -\ 
A.HFTHK)); 

voll = AinI*A.LENGTH+ACF*A.CFTHK; 

voll = AAPPBEDVOID*PI*((A.CFRAD-A.CFTHK)*(A.CFRAD-A.CFTHK)-(A.CFRAD- \ 
A.CFTHK-A.PBEDTHK)*(A.CFRAD-A.CFTHK-A.PBEDTHK))*A.LENGTH; 

voll = AoutlI*A.LENGTH + AHF*A.HFTHK; 

transtime=0; 

N = 0; 

for (time = 0; time <= D.t_DELAY+D.delta_t; time += D.delta_t) 


powerdensity = D.INIT_PRDEN; 

W =WI1; 

Qbed = powerdensity * 1000000000 * PI * ((A.CFRAD - A.CFTHK) * \ 
(A.CFRAD - A.CFTHK) - (A.CFRAD - A.CFTHK - A.PBEDTHK) * \ 
(A.CFRAD - A.CFTHK - A.PBEDTHK)) * A.LENGTH; 

Hinl1 = h2H(D.INIT_TIN); 

Pinl = D.INIT_PIN; 

Pin = D.INIT_PIN; 

PinlI = D.IINIT_POUT; 

Pout = D.INIT_POUT; 

Tin = D.INIT_TIN; 


Tout = Tin; 

do 

{ 
W=WI; 
Hout! = Qbed/W + Hinl1; 
do 


{ 
Htemp = h2H(Tout); 
Herror = Hout! - Htemp; 
Tout = Tout + Herror/100000; 
} while ( Herror/10000 > .001 ); 
/* Set up all variables */ 
TbarlI = .5 * (Tin + Tout); 
Pbarl = .5 * (PinI + PinlI); 
PbarlI = .5 * (PinII + PinlII); 
PbarlIII = .5 * (PinII + Pout); 
pbarl = h2density(Tin, PbarI); 
pbarllI = h2density(TbarlII, PbarIT); 
pbarllI= h2density(Tout, PbarlIII); 
uin = h2viscosity(Tin); 
ubarlI = h2viscosity(TbarlI); 
uout = h2viscosity(Tout); 
pinI = h2density(Tin, Pinl); 
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pinll = h2density(Tin, PinI); 
pinllI = h2density(Tout, PinIII); 
pout = h2density(Tout, Pout); 
/* Inlet Resistance Calculations */ 
Rel = W*Deq_in / (uin* Ain); 
RIf = .138*pow(Rel, -.151)*(A.LENGTH/Deq_in) / (2*pbarI*AinI* Ain); 
Ria = (1/(Ainl*AinII) - 1/(AinI*AinI)) / (2*pbar!); 
RIk = (1-(AinI/AoutI))*(1-(AinI/Aoutl)) / (2*pbarI*AinI*AinI) \ 
+ .5 * (1-(ACF/AoutI))*(1-(ACF/Aoutl)) / (2*pinII*ACF*ACF) \ 
+ .5 * (1-(AinII/ACF))*(1-(AinII/ACF)) / (2*pinII* AinII* Ain); 
RIm = A.MFIN / (pbarI* AinI* Ain); 
RIm = RIm + A.MFIN/(pbarl*AoutI* Aoutl); 
RCF = 150*uin*(1-A.CFPOR)*(1-A.CFPOR)/(A.CFDp*A.CFDp*pinlI* Aoutl* W* \ 
A.CFPOR*A.CFPOR*A.CFPOR) \ 
+ 1.75*(1-A.CFPOR)/(A.CFDp*pinII* AoutI*AoutI*A.CFPOR*A.CFPOR* \ 
A.CFPOR); 
RCF = RCF * A.CFTHK; 
RI = RIf/2 + Rla + Rik + RIm + RCF; 
/* Resistance Calculations for the Particle Bed af | 
RI = 150*ubarll*(1-A.PBEDVOID)*(1-A.PBEDVOID)/(A.PBEDDp*A.PBEDDp* \ 
pbarlI*AbarlI*W*A.PBEDVOID*A.PBEDVOID) \ 
+ 1.75*(1-A.PBEDVOID)/(A.PBEDDp*pbarlI*A barII*Abarll*¥A.PBEDVOID); 
RII = RI * A.PBEDTHK; 
Rila = ( 1/(pinIII*AinIT1) - 1/(pinII*AinII) )*((Ain{I+AinII1)/ \ 
(2* AinII* AinlIII)); 
RII = RII + Rifa; 
/* Resistance for Exit Plenum wal | 
Relll = W*Deq_out / (uout* AoutlII); 
RIIIf = .138*pow(Relll, -.151)*(A.LENGTH/Deq_out)/(2*pbarlIII* AoutlII \ 
* Aoutilll); 
Rilla = (1/(AoutIN*AoutHI)-1/(AoutII* Aoutl))/(2*pbarlIy); 
RIIk = .5*(1-AHF/AoutII)*(1-AHF/AoutlII)/(2*pinII]*AHF*AHF) \ 
+ (1-AHF/AinIi1)*(1-AHF/AinII1)/(2*pinIlI*AHF*AHF) \ 
+ (1-AoutITI/AinHI)*(1-AoutI/Aini)/(2*pbarlil* \ 
AoutII* AoutlIIT); 
mfout = AMMFOUT; 
RIUIm = mfout / (pbarIl/*AoutIII* Aout) + mfout / (pbarHI\ 
* AinII* Ainill); 
RIT = RUOIf/2 + Ria + RIIk + RIim; 
RTOTAL = RI + RII + RII; 
W1 =(PinI - Pout)*1000/RTOTAL; 
W1 =sqrt(W1); 
Werror = W - W1; 
Werror = fabs(Werror); 
PDROPI1=RI* W1*W1/1000; 
PDROPII1=RII*W1*W1/1000; 
PDROPIII1=RUI*W1*W 1/1000; 
PinlI=PinI-PDROPI1; 
PinlJ=Pout+PDROPIII1; 
Pbarl=.5*(PinI+PinlI); 
PbarlI=.5*(PinlI+PinlII); 
PbarlII=.5*(PinIII+Pout); 
vout = W1/(pout*AoutlII); 
power = powerdensity; 
} while(Werror > .0000001); 
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HinI2=Hinl1; 
Hinll1=Hin!l1; 
Hinli2=Hinlt1; 
HinliI1=Hout1; 
Hinlli2=Hinlit1; 
Hout2=Houtl; 
powerout = W*(Hout2-Hinl2); 
DD = fmod(N, EE); 
(op == 0) 
{ 
W1=W1* 1000; 
fprintf(out,"%.4NZNANZANANANANANDNANADNIZIn", \ 
time, power, W1, Tin, TbarII, Tout, powerout, \ 
Pinl, PinlII, PinlJI, Pout, vout); 
printf(“time=%NTin=ZNTout=ZNW=%fn", time, Tin, Tout, W1); 
W1=W1/1000; 


} 
DD =: 
N += 1.0; 
if (Tout > 2700) 
goto meltdown; 


Pinl1 = Pin; 

PinlII1 = Pinll; 

PinlI1 = PinlII; 

Pout1 = Pout; 

PbarI1 = Pbarl; 

PbarlI1 = PbarIl; 

PbarlII1 = PbarII; 

mal = ((C.L3RAD*C.L3RAD*#C.L3RAD)-(C.L2RAD*C.L2RAD*C.L2RAD))*C.L3DEN/\ 
(3*C.L3RAD*C.L3RAD); 

ma2 = ((C.L2RAD*C.L2RAD*#C.L2RAD)-(C.LIRAD*C.LIRAD*C.LIRAD))*C.L2DEN/\ 
(3*C.L3RAD*C.L3RAD); 

ma3 = ((C.LIRAD*C.LIRAD*C.LIRAD)- 

(C.FUELRAD*C.FUELRAD*C.FUELRAD))*C.LIDEN/\ 

(3*C.L3RAD*C.L3RAD); 

ma4 = (C.FUELRAD*C.FUELRAD*C.FUELRAD)*C.FUELDEN / (3*C.L3RAD*C.L3RAD); 

mat = mal+ma2+ma3+ma4; 

pfuelave=(1/( 1.3333*PI*pow(C.L3RAD,3) ) )*( C-FUELDEN*pow(C.FUELRAD,3)\ 
+ C.LIDEN*( pow(C.L1RAD,3)-pow(C.FUELRAD,3) ) + C.L2DEN*( pow(C.L2RAD,3)\\ 
-pow(C.LIRAD,3) ) + C.L3DEN*( pow(C.L3RAD,3)-pow(C.L2RAD,3) ) ); 

Cpbar = (mal*C.L3Cp + ma2*C.L2Cp + ma3*C.L1Cp + ma4*C.FUELCp)/mat; 

U1 = C.L3RAD*C.L2RAD*C.L3K/(C.L3RAD*C.L3RAD*(C.L3RAD-C.L2RAD)); 

U2 = C.L2RAD*C.LIRAD*C.L2K/(C.L3RAD*C.L3RAD*#(C.L2RAD-C.LIRAD)); 

U3 = C.LIRAD*C.FUELRAD*C.LIK/(C.L3RAD*C.L3RAD*(C.LIRAD-C.FUELRAD)); 

U4 = 2*C.FUELK/(3*C.L3RAD); 

UT = 1/1/U1 + 1/02 + 1/03 + 1/04); 

fea T/T: 

f2 = UT/U2 + f1; 

f3 = UT/U3 + f2; 

fbar = (mal*C.L3Cp*f1 + ma2*C.L2Cp*(f1+f2) + ma3*C.L1Cp*(f2+f3) +\ 

ma4*C.FUELCp*(f3+1))/(2*mat*Cpbar); 

Vcv=volll/A.PBEDV OID; 

Av=6*(1-A.PBEDVOID)/A.PBEDDp; 

TinJ=Tin; 

Tfl = Tbarll; 

TinJ{=Tinl; 
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TinlII=Tout; 

deltat=D.delta_t; 

i= OO: 

for (time=D.t_DELAY; time < (D.t_DELAY+D.DUR_PRDEN+D.t_AFTER); \ 
time += deltat) 

{ 


powerdensity = D.INIT_PRDEN + ((D.FIN_PRDEN - D.INIT_PRDEN) \ 
* transtime /D.DUR_PRDEN); 
if (powerdensity >= D.FIN_PRDEN) 
powerdensity = D.FIN_PRDEN; 
k = h2heatxfer(TbarI]); 
Pr = h2prandle(TbarII); 
Rell = (W*6)/(AbarlI*A.PBED VOID *ubarll*Av); 
Nu = .8*pow(Rell, .7)*pow(Pr, .33); 
h = Nu*k*2/A.PBEDTHK; 
UBAR = UT*h/(fbar*h+UT); 
power = powerdensity; 
Qs = power* 1000000000/A v; 
BB=mat*Cpbar* Vcv*Av; 
AA=UBAR* Vcv* Av; 
if (transtime == Q); 
Tfl=Qs*Vcv*Av/(AA) + Tharl]; 
if (transtime > 1) 


if (D.DUR_TIN > 0) 


{ 
Tinl=D.INIT_TIN+(D.FIN_TIN - D.INIT_TIN)*(transtime/D.DUR_TIN); 
Hinl2=h2H(Tinl); 


} 
if (transtime >= D.DUR_TIN) 


{ 
TinI = D.FIN_TIN; 
Hinl2 = h2H(Tinl); 
} 


} 
if (D.DUR_PIN > 0) 
PinI1=D.INIT_PIN+(D.FIN_PIN-D.INIT_PIN)*(transtime/D.DUR_PIN); 
if (transtime >= D.DUR_PIN) 
Pinl1] = D.FIN_PIN; 
if (D.DUR_POUT > 0) 
Poutl=D.INIT_POUT+(D.FIN_POUT-D.INIT_POUT)*(transtime/D.DUR_POUT); 
if (transime >= D.DUR_POUT) 
Poutl = D.FIN_POUT; 
/* Control Volume I */ 
M1 = pinl* voll; 
Ii = .5*A.LENGTH/AinI+ (A.OUTRAD-A.CFRAD)/Aout] + A.CFTHK/ACF; 
W2=(W*I1+deltat*(PinI1-PinII1)*1000)/(11+deltat*RI* W); 
HinlI2 = (HinII1 + (deltat/M1)*((volf*(Pbarl1-PbarI)* 1000/deltat) +\ 
W*Hinl2 + 0.00*Qbed))/(1+deltat* W/M1); 
/* Control Volume II */ 
Tf2 =(Tfl + deltat*(Qs*Vcv*Av + AA*Tbarll)/BB)/(1 + deltat*AA/BB); 
Qbed = AA*(Tf2-TbarlII); 
Til =: 
M2 = (pbarlII*vollI + pfuelave*volli*(1-A.PBEDVOID))/2; 
I2 = A.PBEDTHK/AbarI; 
W2=(W*12+deltat*(PinII1-PinlII1)*1000)/(12+deltat* RII*W); 
HinllI2 = (HinII11 + (deltat/(M2)*((Vcv*(PbarII1 -PbarII)* 1000/deltat) \ 
+ W*Hinll2 + 1.00*Qbed ))/(1+deltat*W/M2); 
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/* Control Volume IU */ 
M3 = pinllI*vollll; 
I3 = A.HFTHK/AoutlI + (A.CFRAD-A.CFTHK-A.PBEDTHK \ 
-A.HFTHK)/AinllI + .5*A. LENGTH/Aouull; 
W2=(W*13+deltat*(PinlII1-Pout1)* 1000)/(13+deltat*RIII* W); 
Hout2 = (Hout1 + (deltat/M3)*( (volIII*(PbarIII1-PbarIII)* 1000/deltat) \ 
+ W*HinllI2 + 0.00*Qbed ))/(1+deltat* W/M3); 
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Pbarl = PbarI1; 
PbarlII = PbarII1; 
PbarllI = PbharIII1; 
PinI = Pin!1; 

Pout = Pout1; 
PinU = PinIi1; 
Pin = PinflI1; 


Htemp = h2H(Tinl); 

Herror = Hinl2 - Htemp; 

TinI = TinI + Herror/100000; 
} while ( Herror/100000 > .001 ); 


do 
{ 
Htemp = h2H(Tout); 
Herror = Hout2 - Htemp; 
Tout = Tout + Herror/100000; 
} while ( Herror/100000 > .001 ); 
do 


Htemp = h2H(Tinl]); 
Herror = HinlI2 - Htemp; 
Tinll = TinlI + Herror/100000; 
} while ( Herror/100000 > .001 ); 
do 


Htemp = h2H(TinII)); 

Herror = HinlIQ2 - Htemp; 

TinlII = Tin + Herror/100000; 
} while ( Herror/100000 > .001 ); 


Tbarl=.5*(TinI+TinlI); 
Tbarll=.5*(TinII+TinlII); 
TbarlI0=.5*(TinHI+Tout); 

HinI1 = Hinl2; 

Hinll1 = HinlI2; 

HinII1 = HinlI12; 

Hout1 = Hout2; 

pinI = h2density(Tinl, Pinl); 

pinll = h2density(Tinl, PinIt); 
pinllI = h2density(TinllI, Pini); 
pout = h2density(Tout, Pout); 
pbarI = h2density(Tbarl,Pbar!); 
pbarII = h2density(TbarlIl, PbarII); 
pbarllI = h2density(Tout, PbarIII); 
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uin = h2viscosity(Tbarl); 

ubarII = h2viscosity(TbarI]); 

uout = h2viscosity(T bar); 

/* Inlet Resistance Calculations | 

Rel = (W*Deq_in)/(uin* Ainl); 

RIf = .138*pow(Rel, -.151)*(A.LENGTH/Deq_in)/(2*pbarl* AinI* Ain!); 

Ria = (1/(AinII* Ainld) - 1/(AinI* Ainl))/(2*pbarl); 

RIk = (1-(AinI/Aoutl))*(1-(AinI/Aoutl)) / (2*pbarl*AinI*AinI) = \ 
+ .5 * (1-(ACF/Aoutl))*(1-(ACF/AoutI)) / (2*pinII*ACF*ACF) \ 
+ .5 * (1-(AinII/ACF))*(1-(AinII/ACF)) / (2*pinll* Aint* Aint); 

RIm = A.MFIN / (pbarI* AinI*AinlI); 

RIm = RIm + A.MFIN/(pbarI* AoutI* Aout); 

RCF = 150*uin*(1-A.CFPOR)*(1-A.CFPOR)/(A.CFDp*A.CFDp* pinII* Aoutl* W* \ 
A.CFPOR*A.CFPOR*A.CFPOR) \ 
+ 1.75*(1-A.CFPOR)/(A.CFDp* pinII* AoutI* AoutI*A.CFPOR* A.CFPOR* \ 
A.CFPOR); 

RCF = RCF * A.CFTHK; 

RI = RIf/2 + Rla + Rik + RIm + RCF; 

/* Resistance Calculations for the Particle Bed i 

RIT = 150*ubarlI*(1-A.PBEDVOID)*(1-A.PBEDVOID)/(A.PBEDDp*A.PBEDDp* \ 
pbarlII* AbarlI* W*A.PBEDVOID*A.PBEDVOID)\ 
+ 1,.75*(1-A.PBEDVOID)/((A.PBEDDp*pbarlII* AbarII*AbarlI* A.PBEDVOID); 

RII = RIT * ALPBEDTHK; 

Rila = ( 1/(pmHI* Ain) - 1/(pintI*AinIt) )*((AintI+AinII)/ \ 


(2* AinII* AinII)); 
RI = RIT + Rila; 
/* Resistance for Exit Plenum na 


RellI = W*Deq_out / (uout* AoutlII); 

RIUIf = .138*pow(Relll, -.151)*(A. LENGTH/Deg_out)/(2* pbarIII* AoutIII \ 
* AoutIII); 

Rilla = (1/(AoutIM* AoutIT)-1/(AoutII* AoutII))/(2* pbarIII); 

RITIk = .5*(1-AHF/AoutII)*(1-AHF/A outID)/((2*pinIN*AHF*AHF) \ 
+ (1-AHF/AiniT)*(1-AHF/AinI1)/(2*pinIlI*AHF*AHF) \ 
+ (1-AoutIII/AinIII)*(1-AoutiI/AiniII)/(2*pbarIIi* \ 
AoutIII*AoutiII); 

mfout = AMFOUT; 

Riiim = mfout / (pbarIII*A outIII* Aout) + mfout / (pbarlIII\ 
* AinIII* AinJIJI); 

RIT = RUAf/2 + RUIa + RUIk + Rim; 

/* TOTALS “27 

RTOTAL = RI + RII + RIII; 

W1 = (PinI - Pout)*1000/RTOTAL; 

W1 = sqrt(W1); 

W=WI1; 

PDROPI1 = RI* W*W/1000; 

PDROPII1 = RIJ*W*W/1000; 

PDROPIII1 = RITI*W*W/1000; 

PinII1 = PinI1 - PDROPI1; 

PinlII1 = Poutl + PDROPITI1; 

Pbarl1 = .6*(PinI1+PinII1); 

PbarlI1 = .§6*(PinII1+PinIIJ1); 

PbharlII1 = .6*(PinIIZ1+Pout1); 

vout = W/(pout* AoutIID); 

powerout = W*(Hout2-Hinl2); 

DD = fmod(N, EE); 

if (DD == 0) 

{ 
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W = W*1000: 
fprintf(out,"%.4N %.4N ZN %.3N O.3N %.3N G.2N %.2N %.2N\ 
J%.2N %.2N WIN %.3fn", \ 
time, power, W, Tinl, Tbarll, Tout, powerout, \ 
Pinl1, PinII1, PinIi1, Poutl, vout, Tf1); 
pnntf(tume=ZNTinI=ZNTout=ZNW=%frn", time, TinI, Tout, W); 
W = W/1000; 


} 
pp =(): 
N += 1.0; 
if (Tout > 2700) 
goto meltdown; 
transtime += deltat; 


meltdown: 
if (Tout > 2700) 
{ printf("Tout > 2700"); 
fprintf(out,""Temp exceeds 2700"); 
} 
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void showdefaults( struct fuelelem *e_ptr ) 


int num; 

float entry; 

int Cc; 

_clearscreen(0); 

printf(’THE FOLLOWING ARE PARAMETERS FOR PPBR FUEL ELEMENTSn\n"); 
pnntf("1. Outer Radius = %f m\n", e_ptr->OUTRAD); 

printf("2. Cold Frit Radius = %f m\n", e_ptr->CFRAD); 

pnntf('3. Cold Frnt Thickness = %f m\n", e_ptr->CFTHK); 

pnntf("4. Dia of Cold Frit Particle = %.8f m\n", e_ptr->CFDp); 

pnntf("'S. Cold Fnit Porosity = %f\n", e_ptr->CFPOR); 


printf("6. Particle Bed Thickness = %f m\n", e_ptr->PBEDTHK); 
pnintf("'7. Fuel Particle Diameter = %f m\n", e_ptr->PBEDDp); 
pnntf("8. Particle Bed Void Fraction = %f\n", e_ptr->PBEDVOID); 
printf("9. Hot Frit Thickness = %f m\n", e_ptr->HFTHK); 
printf('10. Hot Frit Porosity = %f\n", e_ptr->HFPOR); 
printf("'11. Inlet Manifold Factor = %f\n", e_ptr->MFIN); 
printf("12. Exit Manifold Factor = %f\n", e_ptr->MFOUT); 
pnntf("13. Fuel Element Length = %f m\n", e_ptr->LENGTH); 
printf(‘\n\n\nCHANGE A PARAMETER? (1=Y or 0=N)"); 
scanf( ""%1", &c ); 
iec ——() ) 

retum, 
printf(’\nPLEASE ENTER THE PARAMETER NUMBER AND VALUE (eg 12, 1.705)Nn"); 
scanf( "%1, %f", &num, &entry ); 
switch (num) 


case I: 
e_ptr->OUTRAD = entry; 
break; 

case 2: 
e_ptr->CFRAD = entry; 
break; 

case 3: 
e_ptr->CFTHK = entry; 
break; 
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case 4: 
e_ptr->CFDp = entry; 
break; 

Case 3: 
e_ptr->CFPOR = entry; 
break; 

case 6: 
e_ptr->PBEDTHK = entry; 
break; 

Case 7. 
e€_ptr->PBEDDp = entry; 
break; 

case 8: 
e_ptr->PBEDVOID = entry; 
break; 

case 9: 
e_ptr->HFTHK = entry; 
break; 

case 10: 
e_ptr->HFPOR = entry; 
break; 

case 11: 
e_ptr->MFIN = entry; 
break; 

case 12: 
e_ptr->MFOUT = entry; 
break; 

case 13: 
e_ptr->LENGTH = entry; 
break; 

default: 
printf’\nnTRY ANOTHER PARAMETER\n"); 
_clearscreen(0Q); 
break; 

} 


retum,; 
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void fuelelemsum ( struct fuelelem *e_ptr, struct fuelelem *f_ptr ) 


_clearscreen (0); 

printf ("PARAMETER\NNELEMENT A\NBASELINE\n\n"); 

printf("1. Outer RadiusWWN@AN%f m\n", e_ptr->OUTRAD, f_ptr->OUTRAD); 
printf("2. Cold Frit Radius\WN%AN%F m\n", e_ptr->CFRAD, f_ptr->CFRAD); 

printf("3. Cold Fnt Thickness\WN%f\N%f m\n", e_ptr->CFTHK, f_ptr->CFTHK); 
printf("4. Dia of Cold Frit ParticleN%.8NN%.8f m\n", e_ptr->CFDp, f_ptr->CFDp); 
printf("5. Cold Frit Porosity\WN%fN%f\n"," e_ptr->CFPOR, f_ptr->CFPOR); 

printf(6. Particle Bed thicknessSW%fANN%f m\n", e_ptr->PBEDTHK, f_ptr->PBEDTHK); 
printf("7. Fuel Particle Diameter ZINN%f m\n", e_ptr->PBEDDp, f_ptr->PBEDDp); 
printf("8. Particle Bed Void Fraction ZfNN%in", e_ptr->PBEDVOID, f_ptr->PBEDV OID); 
printf("9. Hot Frit Thickness\N%ANI%f m\n", e_ptr->HFTHK, f_ptr->HFTHK); 
printf("10.Hot Frit Porosity\WN%f\NZfn", e_ptr->HFPOR, f_ptr->HFPOR); 
printf("11.Inlet Manifold Factort@fAWN%fn"," e_ptr->MFIN, f_ptr->MFIN); 
printf("12.Exit Manifold FactonNZAN%fn", e_ptr->MFOUT, f_ptr->MFOUT); 
printf("13.Length of Element\WN%f m\N%fin", e_ptr->LENGTH, f_ptr->LENGTH); 
printf('\nnHIT ANY KEY TO CONTINUE"); 

getch(); 
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void showfuelpart ( struct fuelpart *g_ptr ) 


int num; 
float entry; 
int C; 
_clearscreen(0); 
pnntf("THE FOLLOWING ARE PARAMETERS FOR THE FUEL PARTICLE\n"); 
prntf("1. Fuel Material = %s\n", g_ptr->FUELMAT); 
pnntf("2. Layer 1 Material = %s\n", g_ptr->L1MAT); 
printf("3. Layer 2 Material = %s\n", g_ptr->L2MAT); 
printf('4. Layer 3 Material = %s\n", g_ptr->L3MAT); 
printf("5. Radius of Fuel = %f m\n", g_ptr->FUELRAD); 
pnntf("6. Radius of Layer 1 = %f m\n", g_ptr->L1RAD); 
pnintf("7. Radius of Layer 2 = %f m\n", g_ptr->L2RAD); 
printf("8. Radius of Layer 3 = %f m\n", g_ptr->L3RAD); 
printf("9. Density of Fuel = %f kg/m3\n", g_ptr->FUELDEN); 
pnntf("10.Density of Layer 1 = %f kg/m3\n", g_ptr->L1DEN); 
printf("11.Density of Layer 2 = %f kg/m3\n", g_ptr->L2DEN); 
printf("12.Density of Layer 3 = %f kg/m3Nn", g_ptr->L3DEN); 
pnntf("13. Cp of Fuel = %f J/kgK\n", g_ptr->FUELCp); 
printf("14. Cp of Layer 1 = %f J/‘kgK\n", g_ptr->L1Cp); 
prntf("15. Cp of Layer 2 = %f J/kgK\n", g_ptr->L2Cp); 
printf("16. Cp of Layer 3 = %f J/kgK\n", g_ptr->L3Cp); 
printf("17. k for Fuel = %f W/m2K\n", g_ptr->FUELK); 
printf("18. k for Layer 1 = %f W/m2K\n", g_ptr->L1K); 
prntf("19. k for Layer 2 = %f W/m2K\o", g_ptr->L2K); 
pnntf("20. k for Layer 3 = %f W/m2K\n", g_ptr->L3K); 
printf(‘\nCHANGE A PARAMETER? (1= Y or O=N)"); 
scanf( "%1", &c ); 
if(c==0 ) 

retum; 
printf(\nENTER PARAMETER NUMBER ( NOT 1-4 ), VALUE (eg 20, 40)\n"); 
scanf( "%i, %of", &num, &entry ); 
switch (num) 


case 5: 
g_ptr->FUELRAD = entry; 
break; 

case 6: 
g_ptr->LIRAD = entry; 
break; 

case 7: 
g_ptr->L2RAD = entry; 
break; 

case 8: 
g_ptr->L3RAD = entry; 
break; 

case 9: 
g_ptr->FUELDEN = entry; 
break; 

case 10: 
g_ptr->L1DEN = entry; 
break; 

case 11: 
g_ptr->L2DEN = entry; 
break; 
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case 12: 
g_ptr->L3DEN = entry; 
break; 

case 13: 
g_ptr->FUELCp = entry; 
break; 

case 14: 
g_ptr->L1Cp = entry; 
break; 

Case 13: 
g_ptr->L2Cp = entry; 
break; 

case 16: 
g_ptr->L3Cp = entry; 
break; 

case 17: 
g_ptr->FUELK = entry; 
break; 

case 18: 
g_ptr->LIK = entry; 
break; 

case 19: 
g_ptr->L2K = entry; 
break; 

case 20: 
g_ptr->L3K = entry; 
break; 

default: 
prntf("CANNOT CHANGE THIS PARAMETER...SORRY An"); 
getch(); 
break; 

} 
retum; 
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void showconditions ( struct conditions *h_ptr ) 


{ 


int num; 

float entry; 

int C; 

_clearscreen(0); 

printf("THE FOLLOWING INITIAL AND FINAL CONDITIONS ARE SET\n\n"); 
printf("1. Initial Inlet Temperature = %f K\n", h_ptr->INIT_TIN); 


printf("2. Final Inlet Temperature = %f K\n", h_ptr->FIN_TIN); 

printf("3. Duration of Temperature Change = %f sec\n", h_ptr->DUR_TIN); 
printf(""4. Initial Inlet Pressure = %of KPa\n", h_ptr->INIT_PIN); 

printf('S. Final Inlet Pressure = %f KPa\n", h_ptr->FIN_PIN); 

printf("6. Duration of Pressure Change = %f sec\n", h_ptr->DUR_PIN); 
printf(’'7. Initial Outlet Pressure = %f KPa\n", h_ptr->INIT_POUT); 
printf(’8. Final Outlet Pressure = %f KPa\n", h_ptr->FIN_POUT); 
prntf("9. Duration of Pressure Change = %f sec\n", h_ptr->DUR_POUT), 
printf("10. Initial Power Density = %f GW/m3\n", h_ptr->INIT_PRDEN); 
printf("11. Final Power Density = %f GW/m3\n"," h_ptr->FIN_PRDEN); 


printf("12.Duration of Power Density Change = %f sec\n", h_ptr->DUR_PRDEN), 
printf("13. Time Delay to Transient = %f sec\n", h_ptr->t_DELAY); 
printf("14. Time After Transient = %f sec\n", h_ptr->t_AFTER); 

printf("15. Time Step = %f sec\n", h_ptr->delta_t); 

printf("16. Initial Guess at Mass Flowrate = %f kg/sec\n", h_ptr->INIT_W); 
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printf \nCHANGE A PARAMETER? (1=Y or 0=N)"); 
scanf( "%i", &c ); 
if (c==0 ) 
retum; 
printf(\nENTER PARAMETER NUMBER, VALUE (eg 20, 40)\n"); 
scanf( "%i, Yof", &num, &entry ); 
switch (num) 
{ 
case 1: 
h_ptr->INIT_TIN = entry; 
break; 
case 2: 
h_ptr->FIN_TIN = entry; 
break; 
case 3: 
h_ptr->DUR_TIN = entry; 
break; 
case 4: 
h_ptr->INIT_PIN = entry; 
break; 
case 5: 
h_ptr->FIN_PIN = entry; 
break; 
case 6: 
h_ptr->DUR_PIN = entry; 
break; 
case 7: 
h_ptr->INIT_POUT = entry; 
break; 
case 8: 
h_ptr->FIN_POUT = entry; 
break; 
case 9: 
h_ptr->DUR_POUT = entry; 
break; 
case 10: 
h_ptr->INIT_PRDEN = entry; 
break; 
case 11: 
h_ptr->FIN_PRDEN = entry; 
break; 
case 12: 
h_ptr->DUR_PRDEN = entry; 
break; 
case 13: 
h_ptr->t_DELAY = entry; 
break; 
case 14: 
h_ptr->t_AFTER = entry; 
break; 
case 15: 
h_ptr->delta_t = entry; 
break; 
case 16: 
h_ptr->INIT_W = entry; 
break; 
default: 
printf’ CANNOT CHANGE THIS PARAMETER...SORRY.\n"); 
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getch(); 
break; 


returm; 
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double h2density(double xx, double yy) 


float dens; 
dens = .08185 * 300 / xx * yy / 101.325; 
retum dens; 
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double h2viscosity(double xx) 


float visc; 
visc = .000008963 * pow( xx/300, .6733); 
return viSC; 
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double h2heatxfer(double xx) 
{ 


float tem1; 

int 1; 

float tem2; 

float heat[45] = {0,.0362, .0665, .0981, .1282, .1561, .182, .206, .228 
iy) as 2 PAS PRS ol les ase So fo 3511, 3665, .384, .398,\ 
‘412, 426, .44, .452, .464, .476, .488, .500, .512, \ 
24, 536, 548, 560, .572, .584, .596, .608, .62, \ 
.632, .644, .656, .668, .680, .692, .704, .716, .728 }; 

tem2 = xx / 50.0; 

1 = floor(tem2); 

tem1 = heat[i] + (heat[i+1] - heat[i]) * ((xx - (i * 50)) /50 ); 

returm tem 1; 
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double h2cp(double xx) 
{ 


float temcp; 
i ( xx <= 420 ) 
temcp = 4.1868 * 1000 * (1.5395 + .0150825 * xx - 4.02449e-5 *\ 
xx * xx + 3.63544e-8 * xx * xx * xx); 
else 
temcp = 4.1868 * 1000 * (3.58927 - 5.55096e-4 * xx + 6.94235e-7 * xx\ 
Px - 145 1556-10 * xx * xx * xx): 
return temcp; 
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float h2H (float xx) 
{ 


float temH; 
float temH1; 
temH=0; 
temH 1=0; 

if (xx <= 420) 


{ 
temH=4.1868*1000*(1.5395*xx + .0150825*xx*xx/2 - 4.02449e-5* x 
XX*Xx*xx/3 + 3.63544e-8 *¥xx*xx*xx*xx/4); 
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temH=temH - 4.1868*1000*(1.5395*50 + .0150825*50*50/2 - 4.02449e-5* \ 
50*50*50/3 + 3.63544e-8*50*50*50*50/4); 


else 


{ 

temH1=4.1868* 1000*(3.58927*xx - 5.55096e-4*xx*xx/2 + 6.94235e-7*  \ 
XX *xx*xx/3 - 1.45155e-10*xx*xx*xx *xx/4); 

temH1l=temH1 - 4.1868*1000*(3.58927*420 - 5.55096e-4*420*420/2+ \ 
6.9423 5e-7*420*420*420/3 - 1.45155e-10*420*420*420*420/4); 

temH=temH1 + 4.1868*1000*(1.5395*420+.0150825 *420*420/2-4.02449e-5* \ 
420*420*420/3 + 3.63544e-8*420*420*420*420/4); 

temH=temH - 4.1868*1000*(1.5395*50 + .0150825*50*50/2 - 4.02449e-5* \ 
50*50*50/3 + 3.63544e-8*50*50*50*50/4); 


retum temH; 
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float h2prandle(float xx) 
{ 


float tem1; 

int 1; 

float tem2; 

float pran[23] = {.000, .712, .719, .706, .690, .675, .664, \ 
.659, .664, .676, .686, .703, .715, .733, \ 
.748, .763, .778, .793, .808, .823, .838, \ 
853, .868 }; 

tem2 = xx / 100.0; 

i = floor(tem2); 

tem1l = pran[i] + (pran[i+1] - pran[i]) * ((xx - (a * 100)) / 100 ); 

return tem]; 
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APPENDIX E: 


TABLE OF SYMBOLS AND ACRONYMS 


(Fa) peep 


je 


€q 


cross sectional area perpendicular to flow 


particle surface area per unit volume 

smaller areas involved in expansion/contraction 
average particle specific heat 

effective particle diameter 

particle diameter 

equivalent channel diameter 


void fraction or porosity 


equivalent pressure drop due to acceleration 


spatial acceleration losses in the particle bed 


equivalent resistive pressure drop 

equivalent pressure drop due to friction 

equivalent pressure drop due to the manifold effect 
equivalent pressure drop due particle bed effects 


equivalent loss due to sudden expansion/contraction 


TaN, 





If 


friction factor 


heat transfer coefficient 
specific enthalpy 


control volume inertia 


coolant heat conductivity 


expansion/contraction coefficients 

average length of travel through particle bed 
average length of travel through plenums 
particle mass per unit surface area 


mass within the control volume 


particle bed Nusselt number 
wetted perimeter 


average pressure within the control volume 


pressure drop across control volume due to the resistance 


heat energy 
heat generation per surface area 


interphase heat transfer 


Reynolds Number 


total equivalent resistance 
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S,, particle surface area per unit volume 


rs 


effective fuel particle temperature 


T, bulk temperature of coolant 
T; fuel particle surface temperature 


U effective over all heat transfer coefficient 


LL viscosity 


= 


» superficial velocity 
Voy size of the control volume 


W mass flowrate 


ACRONYMS 
ALMCR Advanced Liquid Metal Cooled Reactor 
BNL Brookhaven National Laboratory 
HTGR High Temperature Gas Cooled Reactor 
IMEO Initial Mass in Earth Orbit 
INEL Idaho National Engineering Laboratory 
LEO Low Earth Orbit 
NERVA Nuclear Engine for Rocket Vehicle Application 
NTR Nuclear Thermal Rocket 
PBR Particle Bed Reactor 


PFNTR Pressure Fed Nuclear Thermal Rocket 
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SNAP Space Nuclear Applications Program 


SNL Sandia National Laboratory 
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